sources (version 1.5)
You can also look at:
#! /usr/bin/ksh
###########################################################################
#
#
# SYNTAX: sources
#
# BRIEFLY: Find sources in detector fields of view
#
# DESCRIPTION: This proceedure finds X-ray sources from the totaled sky
# DESCRIPTION: and exposure maps.
# DESCRIPTION: The proceedure then produced FITS format source lists
# DESCRIPTION: and creates regions files for each source and the
# DESCRIPTION: "source 0" region excluding all sources.
# DESCRIPTION: <p>
# DESCRIPTION: The positions and extraction
# DESCRIPTION: radii for the sources are determined in a four step
# DESCRIPTION: process.
# DESCRIPTION: <ol>
# DESCRIPTION: <li> Smooth the summed sky and exposure maps
# DESCRIPTION: with the vsmooth STOOL. This does a variable size
# DESCRIPTION: boxcar smooth which uses the minimum square box of pixels
# DESCRIPTION: which contains at least 50 events. The same size
# DESCRIPTION: box is applied to the counts image and to the exposure map and
# DESCRIPTION: the results are divided to produce a smoothed, exposure-corrected
# DESCRIPTION: image. Pixels with 0 exposure are assigned a pixel value of 0.
# DESCRIPTION: Note that
# DESCRIPTION: each combination of instrument (SIS or GIS) and
# DESCRIPTION: image resolution (256 or possibly 64 for the GIS), has three
# DESCRIPTION: images corresponding to lo (<2keV), hi (>2keV) and all
# DESCRIPTION: energy bands.
# DESCRIPTION:
# DESCRIPTION: <li> Detect sources in the smoothed images.
# DESCRIPTION: This is done using the ascasource STOOL which reports the
# DESCRIPTION: center pixel of each source s well as a suggested extraction
# DESCRIPTION: radius and a maximum extraction radius to avoid confusion
# DESCRIPTION: with other sources. Sources are detected with the following
# DESCRIPTION: algorithm:
# DESCRIPTION: <ol>
# DESCRIPTION: <li> Locate the 20 brightest local peaks in the image.
# DESCRIPTION: A local peak is a pixel which is as bright or brighter
# DESCRIPTION: than the eight pixels surrounding it.
# DESCRIPTION: <li> Examine the peaks to see if they are distinct
# DESCRIPTION: and significant.
# DESCRIPTION: Choose the brightest peak. For each of the
# DESCRIPTION: other peaks, compare the minimum intensity along the
# DESCRIPTION: line segment connecting the peaks
# DESCRIPTION: to the intensity of the lower peak.
# DESCRIPTION: If the lower peak is less than 2.5 times brighter
# DESCRIPTION: than the trough separating it from the brightest peak,
# DESCRIPTION: delete the lower peak.
# DESCRIPTION: repeat this for the second brightest surviving peak,
# DESCRIPTION: etc.
# DESCRIPTION: <li> Determine the maximum extraction radii.
# DESCRIPTION: For each source calculate the distance to the
# DESCRIPTION: the lowest point along a straight line to each other
# DESCRIPTION: source. Set the maximum radius to the nearest
# DESCRIPTION: such trough.
# DESCRIPTION:
# DESCRIPTION: <li> Determine the radial distribution of azimuthally
# DESCRIPTION: averaged intensity around each source.
# DESCRIPTION: Pixels with zero intensity are not included in the
# DESCRIPTION: average. The lowest intensity in the radial
# DESCRIPTION: distrubution is taken as the background flux.
# DESCRIPTION: <li> Delete all sources whose peak intensity is less than
# DESCRIPTION: four times the background intensity.
# DESCRIPTION:
# DESCRIPTION: <li> Determine the "suggested" extraction radius.
# DESCRIPTION: This is set to the minimum radius at which the
# DESCRIPTION: azimuthally averaged intensity is more than 1/5
# DESCRIPTION: of the peak intensity and more than twice the
# DESCRIPTION: background intensity. The suggested radius is
# DESCRIPTION: never greater than the maximum radius.
# DESCRIPTION: </ol>
# DESCRIPTION: <li> Determine the extraction radius for each source.
# DESCRIPTION: The default extraction radius is
# DESCRIPTION: 6 arcmin for GIS and 4 arcmin for SIS.
# DESCRIPTION: If the suggested radius is less than the default radius
# DESCRIPTION: (true for point sources) than the default radius is used.
# DESCRIPTION: Otherwise use the suggested radius. In no case is the
# DESCRIPTION: extraction radius larger than the maximum radius.
# DESCRIPTION:
# DESCRIPTION:
# DESCRIPTION: <li> Reconcile the sources detected in the three energy bands.
# DESCRIPTION: and eliminate duplicates.
# DESCRIPTION: This is done by giving a list of all source positions and radii
# DESCRIPTION: to the sift_sources STOOL. The sources are listed in order of
# DESCRIPTION: decreasing intensity with all of the "all" energy bad sources
# DESCRIPTION: listed first, follwed by the "hi" energy sources, and the
# DESCRIPTION: "lo" energy sources. Sift_sources does the following:
# DESCRIPTION: <ol>
# DESCRIPTION: <li> Take the first source in the list.
# DESCRIPTION: <li> Examine all the other sources in the list.
# DESCRIPTION: <li> If the second source is within the extraction radius
# DESCRIPTION: of the first source, eliminate it.
# DESCRIPTION: <li> If the extraction radii of the two sources overlap,
# DESCRIPTION: reduce the radius of the second source so that
# DESCRIPTION: they just touch.
# DESCRIPTION: <li> Repeat these steps for each of the other remaining
# DESCRIPTION: sources
# DESCRIPTION: </ol>
# DESCRIPTION: </ol>
#
# VERSION: 1.5
#
# HISTORY: 0.0 -> 1.0 10/4/96
# HISTORY: Added code to create a FITS BINTABLE list of sources.
# HISTORY:
# HISTORY: 1.0 -> 1.1 12/4/96
# HISTORY: Now copy all relevant keywords from the totsky image headers
# HISTORY: to the smoothed headers plus adding the SEQNUM keyword.
# HISTORY: This is so that the colorpic tool can read some of this info
# HISTORY: to lable the plot it creates.
# HISTORY: Modified to accomodate the new version of colorpic which
# HISTORY: creates a postscript plot.
# HISTORY: Also stopped deleting the smoothed images.
# HISTORY: Modified the source catalog format by deleting the spectral
# HISTORY: fitting columns and adding the sequence column
# HISTORY:
# HISTORY: 1.1 -> 1.2 12/19/96
# HISTORY: Now generate region files for each source in this routine
# HISTORY:
# HISTORY: 1.2 -> 1.3 1/24/97
# HISTORY: Supplied a missing "let" in calculating the default extraction
# HISTORY: radius for gis data. This resulted in gis sources with the
# HISTORY: default radius being dropped. Also deleted a few unused reads from
# HISTORY: the parfile.
# HISTORY: 1/28/97 supplied code for the background region file in the case
# HISTORY: where there are no foreground sources. Also added comments to the
# HISTORY: tops of each region file.
# HISTORY: 2/4/97 modified the source catalog format.
# HISTORY: 2/7/97 Began adding the same keywords to the sourcecat as to the
# HISTORY: last smoothed image
# HISTORY:
# HISTORY: 1.3 -> 1.4 3/26/97
# HISTORY: began clipping smoothed image pixels with low but non-zero
# HISTORY: exposures. This was to eliminate the spurrious one-pixel
# HISTORY: sources which can crop up around the edges of GIS images
# HISTORY: of very bright sources.
# HISTORY: Note, that this requires STOOLS.3.2 or higher. Also ascasource
# HISTORY: in STOOLS.3.2 corrects a pixel numbering problem in earlier versions.
# HISTORY: pixels are now numbered from "1" instead of "0" in the source
# HISTORY: catalog.
# HISTORY:
# HISTORY: 1.4 -> 1.5 3/31/97
# HISTORY: Now included extended/compact flag in sourcecat. Note this requires
# HISTORY: sift-sources in STOOLS 3.3 or higher.
# HISTORY: changed grouping when calculating defaultr
# HISTORY: Now create region files in detector coordinates
# HISTORY: Source-free region files with no sources are now a BOX region
# HISTORY: instead of a CIRCLE
#
##############################################################################
#DEBUG=1
columns=columns.tmp
data=data.tmp
header=header.tmp
$UTIL/milestone "Detecting sources in summed images"
sourcelist=source.tmp
expoclipfrac=$($UTIL/read_parfile $PARFILE expoclipfrac )
for inst in gis sis; do
########################################
# imaging mode for each instrument
########################################
case "$inst" in
sis) mode=02;;
gis) mode=70;;
*) $UTIL/exception $0 1 "Unknown instrument $inst"
esac
####################
# debugging output
#####################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: inst=$inst
echo ${0##*/}: mode=$mode
fi
##################################
# loop over all image dimensions
##################################
list=$($UTIL/generate_filename totsky $inst any $mode any all )
list=$(ls $list 2>/dev/null)
for imagefile in $list; do
dimen=$($UTIL/read_fits_keyword ${imagefile}[0] NAXIS1 )
####################
# debugging output
#####################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: imagefile=$imagefile
echo ${0##*/}: dimen=$dimen
fi
#########################
# loop over all energies
#########################
rm -f $sourcelist
touch $sourcelist
for nrg in all hi lo; do
totsky=$($UTIL/generate_filename totsky \
$inst $dimen $mode any $nrg )
totexpo=$($UTIL/generate_filename totexpo $inst $dimen $mode any)
smooth=$($UTIL/generate_filename smooth \
$inst $dimen $mode any $nrg )
#########################################################
# read exposure keyword and set exposure clipping level
#########################################################
exposure=$($UTIL/read_fits_keyword $totexpo[0] EXPOSURE )
clip_expo=$($STOOLS/equals "$exposure * $expoclipfrac" )
####################
# debugging output
#####################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: nrg=$nrg
echo ${0##*/}: totsky=$totsky
echo ${0##*/}: totexpo=$totexpo
echo ${0##*/}: smooth=$smooth
echo ${0##*/}: exposure=$exposure
echo ${0##*/}: clip_expo=$clip_expo
fi
############################
# create a smoothed flux map
############################
$UTIL/log_entry "Smoothing $totsky with $totexpo"
$UTIL/log_entry "Clipping exposures below $clip_expo seconds"
$STOOLS/vsmooth $totsky $totexpo $smooth $clip_expo \
>stdout.log 2>stderr.log
$UTIL/stool_test vsmooth $? $0 2 stdout.log stderr.log
####################################################
# copy image header keywords into smoothed image
####################################################
$UTIL/setup_parfile $FTOOLS/fdump.par infile=$totsky \
outfile=out.tmp \
columns=" " \
rows=" " \
prhead=yes \
prdata=no \
showcol=no \
showunit=no \
showrow=no \
showscale=no
rm -f out.tmp
rm -f header.tmp
$FTOOLS/fdump >stdout.log 2>stderr.log
$UTIL/ftool_test fdump $? $0 2 stdout.log stderr.log
grep -v SIMPLE out.tmp |\
grep -v NAXIS |\
grep -v BITPIX |\
grep -v COMMENT|\
grep -v HISTORY|\
grep -v CONTENT|\
grep -v FILIN |\
grep -v CREATOR|\
grep -v END >$header
rm -f out.tmp
sequence=$($UTIL/read_parfile $JOBPAR sequence)
echo "SEQNUM = $sequence" >>$header
$UTIL/setup_parfile $FTOOLS/fmodhead.par infile=${smooth}[0] \
tmpfil=$header
$FTOOLS/fmodhead >stdout.log 2>stderr.log
$UTIL/ftool_test fmodhead $? $0 2 stdout.log stderr.log
#####################################
# detect sources in smoothed image
####################################
$UTIL/log_entry "Detecting sources in $smooth"
$STOOLS/ascasource $smooth >stdout.log 2>stderr.log
code=$?
cat stdout.log >>$sourcelist
$UTIL/stool_test ascasource $? $0 2 stdout.log stderr.log
done
#############################################################
# default radius based on instrument and image resolution
############################################################
case "$inst" in
sis)
let defaultr=38*$dimen/320
let compactr=10*$dimen/320
;;
gis)
let defaultr=24*$dimen/256
let compactr=10*$dimen/256
;;
esac
$UTIL/log_entry "Determining extraction radii"
#########################################################
# determine extraction radius for each source
#########################################################
rm -f z.tmp
touch z.tmp
nsources=$(wc -l <$sourcelist )
typeset -i source=1
while [ $source -le $nsources ]; do
i0=$( awk ' NR == line {print $1}' line=$source $sourcelist )
j0=$( awk ' NR == line {print $2}' line=$source $sourcelist )
maxr=$(awk ' NR == line {print $4}' line=$source $sourcelist )
minr=$(awk ' NR == line {print $5}' line=$source $sourcelist )
##################
# determine radius
###################
radius=$defaultr
if [ $minr -gt $defaultr ]; then
radius=$minr
fi
if [ $radius -gt $maxr ]; then
radius=$maxr
fi
##################################
# is it extended?
##################################
if [ $minr -gt $compactr ]; then
extended=T
else
extended=F
fi
##################################
# record the position and radius
##################################
echo "$i0 $j0 $radius $extended" >> z.tmp
let source=$source+1
done
rm -f $sourcelist
##################################################
# sift through the sources to eliminate overlaps
##################################################
$UTIL/log_entry "Eliminating redundant sources"
$STOOLS/sift_sources <z.tmp >stdout.log 2>stderr.log
code=$?
cp stdout.log $sourcelist
$UTIL/stool_test sift_sources $code $0 2 stdout.log stderr.log
rm -f z.tmp
#########################
# make a nice picture
#########################
losky=$($UTIL/generate_filename smooth $inst $dimen $mode any lo )
allsky=$($UTIL/generate_filename smooth $inst $dimen $mode any all )
hisky=$($UTIL/generate_filename smooth $inst $dimen $mode any hi )
sourcepic=$($UTIL/generate_filename sourcepic $inst $dimen $mode any)
rm -f $sourcepic
awk '{print $1" "$2" "$3}' $sourcelist |\
$STOOLS/colorpic $sourcepic $losky $allsky $hisky \
>stdout.log 2>stderr.log
$UTIL/stool_test colorpic $? $0 2 stdout.log stderr.log
########################################
# make a FITS format list of all sources
########################################
sourcecat=$($UTIL/generate_filename sourcecat $inst $dimen $mode any)
$UTIL/log_entry "Creating fits catalog of sources: $sourcecat"
rm -f $sourcecat
rm -f $data
rm -f $columns
#######################################
# column info for source catalog files
######################################
rm -f $columns
echo "sequence J" >> $columns
echo "index I" >> $columns
echo "xpixel I pixels" >> $columns
echo "ypixel I pixels" >> $columns
echo "rpixel I pixels" >> $columns
echo "RA E degrees" >> $columns
echo "Dec E degrees" >> $columns
echo "rdeg E degrees" >> $columns
echo "hardcnts J counts" >> $columns
echo "softcnts J counts" >> $columns
echo "extended L " >> $columns
##############################
# read coordinate scale
##############################
totsky=$($UTIL/generate_filename totsky $inst $dimen $mode any all)
totsky_lo=$( $UTIL/generate_filename totsky $inst $dimen $mode any lo)
totsky_hi=$( $UTIL/generate_filename totsky $inst $dimen $mode any hi)
crpix1=$($UTIL/read_fits_keyword ${totsky}[0] CRPIX1 )
crval1=$($UTIL/read_fits_keyword ${totsky}[0] CRVAL1 )
cdelt1=$($UTIL/read_fits_keyword ${totsky}[0] CDELT1 )
crpix2=$($UTIL/read_fits_keyword ${totsky}[0] CRPIX2 )
crval2=$($UTIL/read_fits_keyword ${totsky}[0] CRVAL2 )
cdelt2=$($UTIL/read_fits_keyword ${totsky}[0] CDELT2 )
################################
# list sources to an ascii file
################################
nsources=$(wc -l <$sourcelist )
touch $data
typeset -i source=1
while [ $source -le $nsources ]; do
x=$( awk ' NR == line {print $1}' line=$source $sourcelist )
y=$( awk ' NR == line {print $2}' line=$source $sourcelist )
r=$( awk ' NR == line {print $3}' line=$source $sourcelist )
extended=$( \
awk ' NR == line {print $4}' line=$source $sourcelist )
ra=$( $STOOLS/equals "($x-$crpix1)*$cdelt1+$crval1" )
dec=$($STOOLS/equals "($y-$crpix2)*$cdelt2+$crval2" )
rdeg=$($STOOLS/equals "abs(${r}*${cdelt1})" )
softcnts=$($STOOLS/counts_in_circle $totsky_lo $x $y $r )
hardcnts=$($STOOLS/counts_in_circle $totsky_hi $x $y $r )
echo "$sequence $source $x $y $r $ra $dec $rdeg" \
"$softcnts $hardcnts $extended" >> $data
let source=$source+1
done
####################################################
# clear some unwanted things from the header file
####################################################
rm -f z.tmp
grep -v RPIX $header |\
grep -v RVAL |\
grep -v DELT |\
grep -v DUNIT |\
grep -v CROTA |\
grep -v CTYPE |\
grep -v CHANMIN |\
grep -v CHANMAX |\
grep -v CHANTYPE |\
grep -v TOTCTS |\
grep -v USER |\
grep -v HDUCLAS |\
grep -v EXPOSURE >z.tmp
mv z.tmp $header
######################################
# create the source catalog file
######################################
$UTIL/setup_parfile $FTOOLS/fcreate.par cdfile=$columns \
datafile=$data \
headfile=$header \
outfile=$sourcecat \
tbltype=binary
$FTOOLS/fcreate >stdout.log 2>stderr.log
$UTIL/ftool_test fcreate $? $0 2 stdout.log stderr.log
rm -f $data
rm -f $columns
rm -f $sourcelist
rm -f header.tmp
done #loop over image dimension
done #loop over instrument
######################################################################
# generate region files for all instruments, dimensions and sources
######################################################################
$UTIL/log_entry "Generating region files"
regiondir=$($UTIL/read_parfile $PARFILE regiondir )
#########################
# loop over instruments
#########################
for inst in s0 s1 g2 g3; do
case "$inst" in
s?)
sourceinst=sis
display=SIS${inst#s}
;;
g?)
sourceinst=gis
display=GIS${inst#g}
;;
esac
#################################
# loop over source catalogs
#################################
list=$($UTIL/generate_filename sourcecat $sourceinst any any any)
list=$(ls $list 2>/dev/null)
for sourcecat in $list; do
dimen=$($UTIL/parse_filename index $sourcecat )
mode=$( $UTIL/parse_filename mode $sourcecat )
nsources=$($UTIL/read_fits_keyword ${sourcecat}[1] NAXIS2 )
################################
# image binning factor
################################
case "$inst" in
s?)
bin=4
let defaultr=38*$dimen/320
;;
g?)
bin=1
let defaultr=24*$dimen/256
;;
esac
############################
# debugging output
############################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: sourcecat=$sourcecat
echo ${0##*/}: inst=$inst
echo ${0##*/}: dimen=$dimen
echo ${0##*/}: mode=$mode
echo ${0##*/}: bin=$bin
echo ${0##*/}: defaultr=$defaultr
echo ${0##*/}: nsources=$nsources
fi
######################################
# header for source-free region
######################################
freeregion=$($UTIL/generate_filename sourceregion $inst $dimen \
$mode none 0 )
rm -f $freeregion
echo "# This is the $display, source 0 region file" >$freeregion
echo "# at $dimen pixel resolution." >>$freeregion
echo "# It gives the region NOT occupied by any of the" >>$freeregion
echo "# detected sources. " >>$freeregion
echo "# This region file is in DETECTOR coordinates" >>$freeregion
case "$inst" in
s?)
echo "# Data extracted from this region may be" >>$freeregion
echo "# used cautiously for background subtraction ">>$freeregion
echo "# of SIS data for which vignetting is small." >>$freeregion
echo "#" >>$freeregion
;;
g?)
echo "# Vignetting makes it unwise to use data" >>$freeregion
echo "# extracted from this region for background" >>$freeregion
echo "# subtraction." >>$freeregion
;;
esac
if [ $nsources -eq 0 ]; then
#######################################################
# no foreground sources, so background is entire field
#######################################################
echo "# There are no foreground sources, so this" >>$freeregion
echo "# region occupies the entire field of view" >>$freeregion
echo "#" >>$freeregion
center=$($STOOLS/equals "$dimen/2")
echo " BOX($center,$center,$center,$center)" >>$freeregion
fi
##########################
# loop over sources
##########################
source=1
while [ $source -le $nsources ]; do
##############################
# source region file
##############################
region=$($UTIL/generate_filename sourceregion $inst $dimen \
$mode none $source )
rm -f $region
echo "# This is the $display, source $source region file">$region
echo "# at $dimen pixel resolution" >>$region
echo "# This region file is in DETECTOR coordinates" >>$region
echo "#" >>$region
x=$($UTIL/read_bintable $sourcecat xpixel $source )
y=$($UTIL/read_bintable $sourcecat ypixel $source )
r=$($UTIL/read_bintable $sourcecat rpixel $source )
###################################################
# convert sky coordinates to detector coordinates
###################################################
let x=$x*$bin
let y=$y*$bin
#############################
# first try 2.0 pixel radius
#############################
result=$($UTIL/sky2det $x.0 $y.0 $inst $mode 2.0 )
if [ "$result" = "No photons wthin radius" ]; then
######################################
# try the extraction radius
######################################
$UTIL/log_entry "No photons in 2.0 pixel radius"
result=$($UTIL/sky2det $x.0 $y.0 $inst $mode $r.0 )
if [ "$result" = "No photons wthin radius" ]; then
##############################################
# no photons in full radius, somethings wrong
##############################################
$UTIL/log_entry \
"No photons for inst $inst, dimen $dimen, source $source"
continue
fi
fi
#######################################################
# extract detector coordinates and record region file
######################################################
detx=$(echo "$result" |awk '{print $1}' )
dety=$(echo "$result" |awk '{print $2}' )
detx=$($STOOLS/equals "int(($detx-1)/$bin + 1.5)" )
dety=$($STOOLS/equals "int(($dety-1)/$bin + 1.5)" )
echo " CIRCLE($detx,$dety,$r)" >>$region
##########################
# source-free region file
##########################
if [ $r -lt $defaultr ]; then
r=$defaultr
fi
echo "-CIRCLE($detx,$dety,$r)" >>$freeregion
let source=$source+1
done # loop over sources
done # loop over source catalogs
done # loop over instruments
rm -f include.tmp
rm -f exclude.tmp
exit 0