images (version 1.7)
You can also look at:
#! /usr/bin/ksh
###########################################################################
#
#
# SYNTAX: images
#
# BRIEFLY: Generate images and exposure maps
#
# DESCRIPTION: This routine creates exposure maps and sky images for all
# DESCRIPTION: of the filtered event files.
# DESCRIPTION: The selection region used in filtering accounted for, but
# DESCRIPTION: the ratefile currently set to "none" in ascaexpo for both SIS
# DESCRIPTION: and GIS files.
#
# VERSION: 1.7
#
# HISTORY: 0.0 -> 1.0 9/19/96
# HISTORY: Now delete individual exposure maps.
# HISTORY:
# HISTORY: 1.0 -> 1.1 9/27/96
# HISTORY: Now make all energy counts map with xselect along with the hi and
# HISTORY: low energy images. This is because ascaexpo does not put in all
# HISTORY: the keywords that the extractor does. Also ascaexpo images have
# HISTORY: floating point pixels instead of integer pixels (in case this
# HISTORY: matters).
# HISTORY:
# HISTORY: 1.1 -> 1.2 10/17/96
# HISTORY: don't generate individual exposure maps for BRIGHT2 event files
# HISTORY: any more. This is because they are not used in the summed maps
# HISTORY: and probably won't be kept anyway. Also they are redundant with
# HISTORY: the corresponding BRIGHT mode files. Before this the BRIGHT2
# HISTORY: exposure maps were not being deleted and weren't included in
# HISTORY: the distribution or trend products.
# HISTORY:
# HISTORY: 1.2 -> 1.3 12/04/96
# HISTORY: Corrected TOTCNTS to TOTCTS for the total counts keyword
# HISTORY: in the image files. This bug was creating a new keyword
# HISTORY: instead of updating the value of the existing one.
# HISTORY: Also started updating instrument keyword to
# HISTORY: GIS2,GIS3 or SIS0,SIS1 as appropriate
# HISTORY:
# HISTORY: 1.3 -> 1.4 2/4/97
# HISTORY: GIF format image from ximage is now exposure corrected.
# HISTORY: 2/5/97 Changed from floatadd to equals STOOLS
# HISTORY: added some quotes to add_fits_keyword calls
# HISTORY:
# HISTORY: 1.4 -> 1.5 3/18/97
# HISTORY: Now ignore event files which have ontime less than a tollerance.
# HISTORY: This is to avoid large regions of very short exposure time which
# HISTORY: can screw up exposure corrected images and source detection.
# HISTORY: Note: this is only done for SIS event files.
# HISTORY:
# HISTORY: 1.5 -> 1.6 3/26/97
# HISTORY: Undid the last change. The same thing has been achieved in
# HISTORY: the sources routine by clipping individual low exposure
# HISTORY: pixels when smoothing and exposure correcting.
# HISTORY:
# HISTORY: 1.6 -> 1.7 5/5/97
# HISTORY: Made compatible with GIS region files which depend on the
# HISTORY: spread discriminator.
#
##############################################################################
#DEBUG=1
$UTIL/milestone "Generating images and exposure maps"
################################################
# loop through all filtered event files
################################################
list=$($UTIL/generate_filename event any any any any )
list=$(ls $list 2>/dev/null)
for event in ${list} ; do
############
# file names
############
inst=$( $UTIL/parse_filename inst $event)
index=$( $UTIL/parse_filename index $event)
mode=$( $UTIL/parse_filename mode $event)
bitrate=$($UTIL/parse_filename bitrate $event)
####################################################################
# skip BRIGHT2 mode files since their time is included in BRIGHT
# mode (converted) files
#####################################################################
if [ "$mode" = "12" ]; then
continue
fi
teldef=$($UTIL/generate_filename teldef $inst $index $mode $bitrate )
rate=$( $UTIL/generate_filename rate $inst $index $mode $bitrate )
expo=$( $UTIL/generate_filename expo $inst $index $mode $bitrate )
rm -f $expo
######################
# debugging output
######################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: inst=$inst
echo ${0##*/}: index=$index
echo ${0##*/}: mode=$mode
echo ${0##*/}: bitrate=$bitrate
echo ${0##*/}: event=$event
echo ${0##*/}: teldef=$teldef
echo ${0##*/}: rate=$rate
echo ${0##*/}: expo=$expo
echo ${0##*/}: sky=$sky
fi
################################################################
###############################################################
##
## Exposure maps
##
##############################################################
##############################################################
$UTIL/log_entry "Generating exposure map $expo"
############################################
# extraction region map for exposure maps
############################################
case "$inst" in
g[23])
rawxbins=$($UTIL/read_fits_keyword ${event}[1] RAWXBINS )
spread=$( $UTIL/read_fits_keyword ${event}[0] S_DSCR )
extractregion=$($UTIL/read_parfile $PARFILE \
${inst}region${rawxbins}$spread )
dir=$($UTIL/read_parfile $PARFILE regiondir )
$UTIL/fetch_file $dir $extractregion
regionmap=${extractregion}_regionmap.tmp
if [ ! -a "$regionmap" ]; then
#####################################
# haven't already generated the map
#####################################
$STOOLS/regionmap $regionmap $extractregion $rawxbins \
>stdout.log 2>stderr.log
$UTIL/stool_test regionmap $? $0 2 stdout.log stderr.log
fi
imath="mul"
;;
s[01])
regionmap="none"
imath="none"
;;
*) $UTIL/exception $0 1 "Unknown instrument $inst"
esac
######################
# debugging output
######################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: imath=$imath
echo ${0##*/}: regionmap=$regionmap
fi
############################
# fetch teldef file
############################
caldir=$($UTIL/read_parfile ./$PARFILE caldir )
$UTIL/fetch_file $caldir $teldef
##########################
# run ascaexpo
##########################
$UTIL/log_entry "running ascaexpo for $event"
$UTIL/setup_parfile $FTOOLS/ascaexpo.par \
datafile=$event \
calfile=$teldef \
attitude=default \
defattpath=./ \
ratefile=none \
satfcol="SAT" \
instfile="$regionmap" \
imath="$imath" \
expofile=$expo \
imapfile=none \
skyfile=none \
attstep=5 \
rebin=-1 \
binoff=0.0
$FTOOLS/ascaexpo >stdout.log 2>stderr.log
$UTIL/ftool_test ascaexpo $? $0 2 stdout.log stderr.log
################################################################
###############################################################
##
## Low, High, and All Energy Images
##
##############################################################
##############################################################
$UTIL/log_entry "Generating low and high energy images for $event"
command=command.tmp
#############################################
# determine PI bin for high/low energy cut
#############################################
case "$inst" in
g*)
pibins=$($UTIL/read_fits_keyword ${event}[0] PHA_BINS)
picut=$($UTIL/read_parfile $PARFILE gis_picut_$pibins);;
s*)
pibins=2048
picut=$($UTIL/read_parfile $PARFILE sis_picut);;
esac
losky=$($UTIL/generate_filename sky $inst ${index} $mode $bitrate lo )
hisky=$($UTIL/generate_filename sky $inst ${index} $mode $bitrate hi )
sky=$($UTIL/generate_filename sky $inst ${index} $mode $bitrate all )
rm -f $losky
rm -f $hisky
rm -f $sky
##############################
# setup xselect command file
##############################
rm -f $command
echo "${event%.*}" >>$command
echo " " >>$command
echo "%ED%off" >>$command
echo "set mission ASCA" >>$command
echo "set datadir ." >>$command
echo "read events $event" >>$command
echo "set XYNAME X Y" >>$command
echo "extract image" >>$command
echo "save image outfile=${sky}" >>$command
echo "filter pha_cut 0 $picut" >>$command
echo "extract image" >>$command
echo "save image outfile=${losky}" >>$command
echo "clear pha_cut" >>$command
echo "filter pha_cut $picut $pibins" >>$command
echo "extract image" >>$command
echo "save image outfile=${hisky}" >>$command
echo "exit save_session=no" >>$command
######################
# run xselect
######################
rm -f ${losky}
rm -f ${hisky}
$FTOOLS/xselect @${command} </dev/null > stdout.log 2> stderr.log
$UTIL/xselect_test $? $0 stdout.log stderr.log
rm -f $command
done
rm -f *_regionmap.tmp
########################################
########################################
##
## Combine Images
##
########################################
########################################
datalist=datalist.tmp
listfile=list.tmp
ppm=pgplot.ppm
export XANADU=$($UTIL/read_parfile $PARFILE xanadu )
export DISPLAY=$($UTIL/read_parfile $PARFILE display)
ppm2gif=$($UTIL/read_parfile $PARFILE ppm2gif )
giftool=$($UTIL/read_parfile $PARFILE giftool )
csh $XANADU/tools/initu.csh
datalist=datalist.tmp
listfile=listfile.tmp
for inst_type in sis gis; do
$UTIL/log_entry "Summing $inst_type images"
###########################################################
# select BRIGHT mode for the SIS and PH mode for the GIS
###########################################################
case "$inst_type" in
sis) mode=02;;
gis) mode=70;;
esac
#################################################################
# list the files and dimensions like this: NAXIS1 filename
#############################################################
list=$($UTIL/generate_filename expo ${inst_type%is}? any $mode any )
list=$(ls $list 2>/dev/null )
rm -f dimenlist.tmp
touch dimenlist.tmp
for file in $list ; do
naxis=$($UTIL/read_fits_keyword ${file}[0] NAXIS1 )
inst=$( $UTIL/parse_filename inst $file )
index=$( $UTIL/parse_filename index $file )
bitrate=$($UTIL/parse_filename bitrate $file )
echo "$naxis ${inst}|${index}|${bitrate}" >>dimenlist.tmp
done
#####################################
# loop through each unique dimension
######################################
for dimen in $(awk '{print $1}'<dimenlist.tmp |sort -n |uniq ); do
rm -f $datalist
touch $datalist
awk '$1 == word {print $2}' word=$dimen dimenlist.tmp \
>$datalist 2>/dev/null
##############################
# debugging output
##############################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: dimen=$dimen
fi
############################################
# loop over the types of files to be summed
############################################
for type in expo all lo hi; do
case "$type" in
expo) summed=$($UTIL/generate_filename \
totexpo $inst $dimen $mode any );;
*)
summed=$($UTIL/generate_filename \
totsky $inst $dimen $mode any $type );;
esac
#######################
# debugging output
#######################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: type=$type
echo ${0##*/}: summed=$summed
fi
##########################
# assemble a list of files
# and sum the exposure
###########################
totexposure=0.0
totcounts=0
rm -f $listfile
for thing in $(cat $datalist ); do
###########################
# determine filename
############################
inst=$( echo $thing |awk -F'|' '{print $1 }' )
index=$( echo $thing |awk -F'|' '{print $2 }' )
bitrate=$(echo $thing |awk -F'|' '{print $3 }' )
case "$type" in
expo) file=$($UTIL/generate_filename \
expo $inst $index $mode $bitrate );;
*)
file=$($UTIL/generate_filename \
sky $inst $index $mode $bitrate $type );;
esac
#######################
# debugging output
#######################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: file=$file
fi
#########################
# echo filename to list
#########################
if [ -a "$file" ]; then
echo "$file" >>$listfile
exposure=$($UTIL/read_fits_keyword ${file}[0] EXPOSURE)
totexposure=$($STOOLS/equals $exposure +$totexposure )
if [ "$type" != "expo" ]; then
counts=$($UTIL/read_fits_keyword ${file}[0] \
TOTCTS )
let totcounts=$totcounts+$counts
fi
else
$UTIL/exception $0 2 "$file missing for image sum"
fi
done # loop over files
########################
# add the images
##########################
rm -f $summed
$UTIL/merge_images $summed $listfile
###########################################
# modify INSTRUME and EXPOSURE keywords
###########################################
case $inst_type in
s*) $UTIL/add_fits_keyword ${summed}[0] INSTRUME "SIS0,SIS1" ;;
g*) $UTIL/add_fits_keyword ${summed}[0] INSTRUME "GIS2,GIS3" ;;
esac
$UTIL/add_fits_keyword ${summed}[0] EXPOSURE "$totexposure"
if [ "$type" != "expo" ]; then
$UTIL/add_fits_keyword ${summed}[0] TOTCTS "$totcounts"
fi
#################################################################
# delete low and high energy individual images and exposure maps
#################################################################
case "$type" in
hi|lo|expo)
for file in $( cat $listfile ); do
rm -f $file
done
;;
esac
rm -f $listfile
done #loop over types
rm -f $datalist
############################
# GIF image
############################
totexpo=$($UTIL/generate_filename totexpo $inst_type $dimen \
$mode any )
totsky=$( $UTIL/generate_filename totsky $inst_type $dimen \
$mode any all )
index=${dimen#dimen}
gif=$($UTIL/generate_filename gif $inst_type $dimen $mode any)
rm -f $gif
$UTIL/log_entry "Running XIMAGE to create $gif"
rm -f $command
echo "read/fits $totsky" >>$command
echo "read/exposure $totexpo" >>$command
echo "smooth" >>$command
echo "cpd /ppm" >>$command
echo "disp/correct/noframe" >>$command
echo "grid" >>$command
echo "scale" >>$command
echo "exit" >>$command
$XANADU/sol/bin/ximage @$command </dev/null >stdout.log 2>stderr.log
$UTIL/stool_test ximage $? $0 2 stdout.log stderr.log
########################
# convert to GIF format
########################
$ppm2gif < $ppm > $gif 2>> stderr.log
$UTIL/exit_test $? $0 2 "while converting a PPM file to GIF format"
rm -f $ppm
rm -f $command
##############################
# interlace the GIF
##############################
$giftool -i -B $gif
$UTIL/exit_test $? $0 2 "while interlacing $gif with $giftool"
done # loop over dimen
rm -f dimenlist.tmp
done # loop over inst
exit 0