lightcurves (version 1.2)
You can also look at:
#! /usr/bin/ksh
###########################################################################
#
#
# SYNTAX: spectra
#
# BRIEFLY: Extract light curves.
#
# DESCRIPTION: Extracts light curves for sources.
# DESCRIPTION: This routine extracts binned light curves (FITS format
# DESCRIPTION: rate files) for each source. The time bin size is set to
# DESCRIPTION: give a fixed mean number of counts per bin. The binned light
# DESCRIPTION: curves are extracted with exposure=1 in xselect so that
# DESCRIPTION: bins must fall entirely within GTIs to be included.
# DESCRIPTION: This routine also creates GIS L1 count rate monitor
# DESCRIPTION: light curves using the FTOOL ghkcurve. The GTIs for those
# DESCRIPTION: curves are created by merging the GTIs in all the filtered
# DESCRIPTION: event files for a given instrument.
#
# VERSION: 1.2
#
# HISTORY: 0.0 -> 1.0 2/10/97
# HISTORY: Now remove lightcurves with fewer than 5 bins.
# HISTORY: Fleshed out the description.
# HISTORY: modified ghkcurve portion to use the new gti_merge utility
# HISTORY:
# HISTORY: 1.0 -> 1.1 3/11/97
# HISTORY: Now extract with exposure=0.5 instead of 1.0
# HISTORY: changed the filename at the top of the plot to the lightcurve
# HISTORY: file name.
# HISTORY: Commented out dead time correction since most of the time it
# HISTORY: doesn't work.
# HISTORY:
# HISTORY: 1.1 -> 1.2 4/14/97
# HISTORY: adapted to use detector coordinate region files
# HISTORY: will now make separate light curves for non-imaging modes.
# HISTORY: Now produce L1 HK light curves for both HIGH bitrate and
# HISTORY: HIGH+MEDIUM bitrates.
#
##############################################################################
#DEBUG=1
$UTIL/milestone "Extracting light curves"
############################
# some temporary file names
############################
eventlist=eventlist.tmp
xselcom=z.xco
command=z.xcm
#########################################
# minimum events for keeping a light curve
#########################################
minlcevents=$($UTIL/read_parfile $PARFILE minlcevents )
meanperlcbin=$($UTIL/read_parfile $PARFILE meanperlcbin )
index=000
#########################
# loop over instruments
#########################
for inst in s0 s1 g2 g3; do
case "$inst" in
s?)
modelist="02 03"
sourcemode="02"
full=sis
;;
g?)
modelist="70 71"
sourcemode=70
full=gis
;;
*)
$UTIL/exception $0 1 "Unknown instrument $inst"
;;
esac
#######################
# debugging output
########################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: inst=$inst
echo ${0##*/}: modelist=$modelist
echo ${0##*/}: full=$full
fi
for mode in $modelist; do
#######################
# debugging output
########################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: mode=$mode
fi
#################################
# set binsize to largest timedel
#################################
binsize=0.0
list=$($UTIL/generate_filename event $inst any $mode any )
list=$(ls $list 2>/dev/null )
#######################################################
# if there are no event files, go on to the next mode
#######################################################
if [ -z "$list" ]; then
continue
fi
for event in $list; do
timedel=$($UTIL/read_fits_keyword ${event}[0] TIMEDEL)
dum=$($STOOLS/floatcmp $timedel $binsize)
############################
# save the largest bin size
############################
if [ "$dum" = "gt" ]; then
minbinsize=$timedel
fi
done
#######################
# debugging output
########################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: minbinsize=$minbinsize
fi
##############################
# init for loop over sources
##############################
case $inst in
s?) dimen=320;;
g?) dimen=$($UTIL/read_fits_keyword ${event}[0] RAWXBINS );;
*) $UTIL/exception $0 1 "Unknown instrument $inst";;
esac
sourcecat=$($UTIL/generate_filename sourcecat \
$full $dimen $sourcemode $bitrate )
################################################
# go to the next mode if there is no sourcecat
################################################
if [ ! -f "$sourcecat" ]; then
$UTIL/log_entry "No source cat $sourcecat"
continue
fi
#########################
# loop over sources
#########################
if [ "$mode" = "71" -o "$mode" = "03" ]; then
#################
# MPC mode
#################
nsources=0
else
nsources=$($UTIL/read_fits_keyword ${sourcecat}[1] NAXIS2 )
fi
######################################################
# only create the source zero light curve if there
# are no sources
#####################################################
if [ $nsources -gt 0 ]; then
typeset -i source=1
else
typeset -i source=0
fi
while [ $source -le $nsources ]; do
lightcurve=$($UTIL/generate_filename lightcurve $inst \
$index $mode any $source)
rm -f $lightcurve
region=$($UTIL/generate_filename sourceregion $inst $dimen \
$sourcemode none \
$source)
#####################################################
# extract an event file first in order to count
# the events and determine a good bin size
####################################################
$UTIL/log_entry "Extracting events from region $region "
$UTIL/log_entry "... and files: $list"
##################################
# generate an xselect command file
##################################
rm -f $xselcom
echo "xsel" >$xselcom
echo "%ED%off" >>$xselcom
echo "set mission ASCA" >>$xselcom
echo "set datadir ." >>$xselcom
for event in $list; do
echo "read events ${event}" >>$xselcom
done
if [ "$mode" != "71" -a "$mode" != "03" ]; then
echo "filter region $region" >>$xselcom
fi
echo "extract events" >>$xselcom
echo "save events outfile=events.tmp use_events=yes" >>$xselcom
echo "exit save_session=no" >>$xselcom
##################################
# feed the comand file to xselect
##################################
rm -f events.tmp
$FTOOLS/xselect @${xselcom} </dev/null >stdout.log \
2>stderr.log
$UTIL/xselect_test $? $0 stdout.log stderr.log
rm -f $xselcom
############################################################
# skip to the next source if no event file was extracted
###########################################################
if [ ! -a events.tmp ]; then
let source=$source+1
continue
fi
##########################################################
# if there were less than a minimum number of events,
# delete the event file and skip to the next iteration
########################################################
totcts=$($UTIL/read_fits_keyword events.tmp[1] NAXIS2 )
if [ $totcts -lt $minlcevents ]; then
$UTIL/log_entry \
"skipping $lightcurve since it would have $totcts events"
rm -f events.tmp
let source=$source+1
continue
fi
#####################################################
# determine bin size
# to produce a set mean number of events per bin
#####################################################
ontime=$($UTIL/read_fits_keyword events.tmp[1] ONTIME )
binsize=$($STOOLS/equals "$ontime * $meanperlcbin / $totcts" )
dum=$($STOOLS/floatcmp $binsize $minbinsize )
if [ "$dum" = "lt" ]; then
binsize=$minbinsize
fi
############################
# debugging output
############################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: totcts=$totcts
echo ${0##*/}: ontime=$ontime
echo ${0##*/}: binsize=$binsize
fi
#################################
# extract the light curve
#################################
$UTIL/log_entry "Extracting $lightcurve with binsize $binsize"
rm -f $xselcom
echo "xsel" >$xselcom
echo "set mission ASCA" >>$xselcom
echo "set datadir ." >>$xselcom
echo "read events events.tmp" >>$xselcom
echo "set binsize $binsize" >>$xselcom
echo "extract curve exposure=0.5" >>$xselcom
echo "save curve ${lightcurve}" >>$xselcom
echo "exit save_session=no" >>$xselcom
##################################
# feed the comand file to xselect
##################################
$FTOOLS/xselect @${xselcom} </dev/null >stdout.log \
2>stderr.log
$UTIL/xselect_test $? $0 stdout.log stderr.log
rm -f $xselcom
rm -f events.tmp
#################################################################
# delete lightcurve and skip ahead if it has less than five bins
#################################################################
naxis2=$($UTIL/read_fits_keyword $lightcurve[1] NAXIS2)
if [ $naxis2 -lt 4 ]; then
$UTIL/log_entry \
"Removing $lightcurve since it has only $naxis2 bins"
rm -f $lightcurve
let source=$source+1
continue
fi
##################################
# dead time correction for GIS data
####################################
case "$inst" in
# g?)
hohoho)
mkf=$($UTIL/any_filename filter )
rm -f mkf.tmp
ls $mkf >mkf.tmp 2>/dev/null
$UTIL/log_entry "Correcting for dead time using:"
$UTIL/file_to_log mkf.tmp
$UTIL/setup_parfile $FTOOLS/ldeadtime.par \
lcfile=$lightcurve \
mkffile="@mkf.tmp"
$FTOOLS/ldeadtime >stdout.log 2>stderr.log
$UTIL/ftool_test ldeadtime $? $0 2 stdout.log stderr.log
rm -f mkf.tmp
;;
esac
###################################
# XRONOS timing analysis and plot
###################################
lcplot=$($UTIL/generate_filename lcplot $inst \
$index $mode any $source)
rm -f $lcplot
rm -f xronos.ps
rm -f pgplot.ps
$UTIL/log_entry "Plotting light curve $lcplot"
rm -f $command
echo "cui command" >>$command
echo "data $lightcurve" >>$command
echo "sta" >>$command
echo "cpd /ps" >>$command
echo "lc1/plt" >>$command
echo "LA File $lightcurve" >>$command
echo "hard /ps" >>$command
echo "quit" >>$command
echo "exit" >>$command
xronos <$command >stdout.log 2>stderr.log
code=$?
###########################################################
# The "Welcome to Xronos" line has
# some wierd non-printing characters in it which must be
# removed
###########################################################
grep -v "Welcome to Xronos" stdout.log |\
grep -v '\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-' >filtout.log
rm -f stdout.log
################################################
# this fortran stuff is not an error
##################################################
grep -v "Note: the following IEEE floating-point" stderr.log |\
grep -v "occurred and were never cleared" |\
grep -v "Inexact;" | grep -v "Underflow;" |\
grep -v "Sun's implementation of IEEE arithmetic is" |\
grep -v "the Numerical Computation Guide." >filt.log
$UTIL/ftool_test xronos $code $0 2 filtout.log stderr.log \
filt.log
mv pgplot.ps $lcplot
rm -f xronos*
rm -f cd_xronos.cm
rm -f $command
let source=$source+1
done #loop over sources
done #loop over mode
done #loop over instrument
#########################################################################
# GIS HK lightcurves:
# These are light curves constructed from the G2_L1 and G3_L1 count
# rate monitors as constructed by ghkcurve
#########################################################################
for inst in g2 g3; do
####################
# debugginbg output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: inst=$inst
fi
#########################################################################
# construct a single GTI table from all of the GIS filtered event files
#########################################################################
rm -f list.tmp
rm -f gti.tmp
list=$($UTIL/generate_filename event $inst any any any )
ls $list 2>/dev/null |awk '{print $1"[2]" }' >list.tmp 2>/dev/null
if [ ! -s "list.tmp" ]; then
########################################################
# skip to next instrument if there are no event files
########################################################
rm -f list.tmp
continue
fi
$UTIL/gti_merge gti.tmp list.tmp
rm -f list.tmp
#########################################
# make a light curve for each HK file
#########################################
for hk in $(ls ft*G${inst#g}HK.fits); do
for bitrate in h hm; do
case "$bitrate" in
h) irate="HIGH";;
hm) irate="HI+MED";;
esac
$UTIL/log_entry "Making L1 light curve of $hk with irate=$irate"
curve=${hk%.fits}${bitrate}.lc
rm -f $curve
$UTIL/setup_parfile $FTOOLS/ghkcurve.par infile=$hk \
outfile=$curve \
gtifile=gti.tmp \
irate="$irate"
$FTOOLS/ghkcurve >stdout.log 2>stderr.log
$UTIL/ftool_test ghkcurve $? $0 2 stdout.log stderr.log
done
done
rm -f gti.tmp
done
exit 0