lightcurves (version 2.1)
You can also look at:
#! /usr/bin/ksh
###########################################################################
#
# SYNTAX: lightcurves
#
# BRIEFLY: Extract light curves for each source
#
# DESCRIPTION: This routine extracts binned light curves (FITS format
# DESCRIPTION: rate files) for each source.
# DESCRIPTION: Each light curve combines data from all of the event files
# DESCRIPTION: for a given instrument - even if mode
# DESCRIPTION: (e.g. discriminator setting) changes between these
# DESCRIPTION: files could affect the count rate.
# DESCRIPTION: <P>
# DESCRIPTION: The time bin size is set to give a mean of
# DESCRIPTION: %meanperlcbin counts per bin, but not to be smaller than
# DESCRIPTION: the coarsest time resolution of the parent event files.
# DESCRIPTION: The binned light
# DESCRIPTION: curves are extracted using exposure=0.5 in <TT>xselect</TT>,
# DESCRIPTION: so that a bin is included only if
# DESCRIPTION: at least half of it falls within a GTI.
# DESCRIPTION: Light curves are only extracted if they contain
# DESCRIPTION: at least %minlcevents events.
# DESCRIPTION: <P>
# DESCRIPTION: The GIS lightcurves are currently <STRONG>NOT</STRONG>
# DESCRIPTION: corrected for dead time due to problems with the current
# DESCRIPTION: (FTOOLS 4.0 and earlier) version of <TT>ldeadtime</TT>.
# DESCRIPTION: <P>
# DESCRIPTION: This routine also creates GIS L1 count rate monitor
# DESCRIPTION: light curves using the FTOOL <TT>ghkcurve</TT>.
# DESCRIPTION: The GTIs for those
# DESCRIPTION: curves are created by merging the GTIs in all the filtered
# DESCRIPTION: event files for a given instrument.
#
# VERSION: 2.1
#
# HISTORY: 0.0 -> 1.0 2/10/97
# HISTORY: Now removes 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 extracts 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 deadtime 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 produces L1 HK light curves for both HIGH bitrate and
# HISTORY: HIGH+MEDIUM bitrates.
# HISTORY:
# HISTORY: 1.2 -> 1.3 6/3/97
# HISTORY: Added check for whether lightcurve was actually extracted.
# HISTORY: Added some additional logging.
# HISTORY: Fixed a bug where minbinsize was not set if there was only one event
# HISTORY: file.
# HISTORY:
# HISTORY: 1.3 -> 1.4 6/9/97
# HISTORY: Added missing space before "]" in test statement.
# HISTORY:
# HISTORY: 1.4 -> 1.5 7/11/97
# HISTORY: Fixed a bug that caused minbinsize to not be calculated correctly.
# HISTORY: Also, changed the logging in the part of the code which calculates
# HISTORY: minbinsize.
# HISTORY:
# HISTORY: 1.5 -> 1.6 8/6/97
# HISTORY: Added a check to skip light curves if the source region filter
# HISTORY: does not exist.
# HISTORY:
# HISTORY: 1.6 -> 1.7 10/9/97
# HISTORY: Now only trys to make ghkcurve light curves for HK files
# HISTORY: for which screened event files have been extracted from the
# HISTORY: corresponding telem file.
# HISTORY:
# HISTORY: 1.7 -> 1.8 10/15/97
# HISTORY: Now deletes light curves which have no counts.
# HISTORY:
# HISTORY: 1.8 -> 1.9 2/27/94
# HISTORY: Stopped deleting light curves with TOTCTS=0, since this is true
# HISTORY: for all light curves.
# HISTORY:
# HISTORY: 1.9 -> 2.0 1998-09-17
# HISTORY: modified to use xanadu initialization
# HISTORY:
# HISTORY: 2.0 -> 2.1 2000-02-29
# HISTORY: Modified to use the lcurve FTOOL to plot light curves
# HISTORY: instead of xronos, which no longer exists and a single
# HISTORY: command driven entity.
# HISTORY: Now explicitly specify detector coordinates for xselect.
#
###########################################################################
#DEBUG=1
$UTIL/milestone "Extracting light curves"
############################
# Some temporary file names
############################
eventlist=eventlist.tmp
xselcom=z.xco
command=lightcurve.pco
###########################################
# 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
#################################
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
minbinsize=0.0
for event in $list; do
timedel=$($UTIL/read_fits_keyword ${event}[0] TIMEDEL)
$UTIL/log_entry "TIMEDEL=$timedel for $event"
############################
# Save the largest bin size
############################
dum=$($STOOL/floatcmp $timedel $minbinsize)
if [ "$dum" = "gt" ]; then
minbinsize=$timedel
fi
done
$UTIL/log_entry "Minimum bin size is $minbinsize seconds"
#############################
# 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)
##############################
# Check if region file exists
##############################
if [ ! -a "$region" ]; then
$UTIL/log_entry \
"Skipping $lightcurve since $region does not exist"
let source=$source+1
continue
fi
################################################
# 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 "set image det" >>$xselcom
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 command file to xselect
###################################
rm -f events.tmp
$FTOOL/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=$($STOOL/equals "$ontime * $meanperlcbin / $totcts")
dum=$($STOOL/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
##################################
$FTOOL/xselect @${xselcom} </dev/null >stdout.log \
2>stderr.log
$UTIL/xselect_test $? $0 stdout.log stderr.log
rm -f $xselcom
rm -f events.tmp
#####################################
# Check if light curve was extracted
#####################################
if [ ! -s "$lightcurve" ]; then
############################
# Light curve not extracted
############################
$UTIL/log_entry "$lightcurve not extracted"
let source=$source+1
continue
fi
#################################################################
# 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 $FTOOL/ldeadtime.par \
lcfile=$lightcurve \
mkffile="@mkf.tmp"
$FTOOL/ldeadtime >stdout.log 2>stderr.log
$UTIL/ftool_test ldeadtime $? $0 2 stdout.log stderr.log
rm -f mkf.tmp
;;
esac
#############################################
# plot the lightcurve using the lcurve FTOOL
#############################################
lcplot=$($UTIL/generate_filename lcplot $inst \
$index $mode any $source)
rm -f $command
echo "dev $lcplot/ps" >>$command
echo "LA File $lightcurve" >>$command
echo "plot" >>$command
echo "quit" >>$command
$UTIL/setup_parfile $FTOOL/lcurve.par nser=1 \
cfile1=$lightcurve \
window="-" \
dtnb=INDEF \
nbint=INDEF \
outfile=" " \
plot=yes \
plotdev=/null \
plotfile="-"
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: plotting light curve
fi
echo "@$command" |\
$FTOOL/lcurve >stdout.log 2>stderr.log
$UTIL/ftool_test lcurve $? $0 2 stdout.log stderr.log
rm -f $command
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: done plotting light curve
fi
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
###################
# Debugging 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
#############################################################
# Extract a list of telemetry files for which there are GTIs
#############################################################
rm -f telem.tmp
for file in $(sed -e's/\[2\]//g' list.tmp); do
telem=$($UTIL/read_fits_keyword $file[0] TLM_FILE)
echo $telem >>telem.tmp
done
rm -f list.tmp
######################################
# Make a light curve for each HK file
######################################
for telem in $(sort telem.tmp | uniq); do
hk=$(echo $telem | sed -e's/\./_/g')
hk=${hk}G${inst#g}HK.fits
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: telem=$telem
echo ${0##*/}: hk=$hk
fi
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 $FTOOL/ghkcurve.par infile=$hk \
outfile=$curve \
gtifile=gti.tmp \
irate="$irate"
$FTOOL/ghkcurve >stdout.log 2>stderr.log
$UTIL/ftool_test ghkcurve $? $0 2 stdout.log stderr.log
done
done
rm -f telem.tmp
rm -f gti.tmp
done
exit 0