spectra (version 2.1)
You can also look at:
#! /usr/bin/ksh
###########################################################################
#
#
# SYNTAX: spectra
#
# BRIEFLY: Extract PI spectra of source and background
#
# DESCRIPTION: This routine extracts spectra for sources, generates response
# DESCRIPTION: matricies, extracts blank sky background spectra, and
# DESCRIPTION: plots the spectra.
# DESCRIPTION: <p>
# DESCRIPTION: Spectra are extracted in sky coordinates from the filtered
# DESCRIPTION: event files. Event files are combined if the spectra extracted
# DESCRIPTION: from them would have the same response matrix and if they can be
# DESCRIPTION: extracted using the same region filter.
# DESCRIPTION: Specifically, SIS data with different CCD modes (1, 2 or 4)
# DESCRIPTION: are kept separate, and GIS data with different numbers of
# DESCRIPTION: PHA bins are kept separate.
# DESCRIPTION: Also, data from different instruments or in different
# DESCRIPTION: modes (BRIGHT BRIGHT2 FAST PH MPC), and data with different
# DESCRIPTION: image resolutions are kept separate.
# DESCRIPTION: The spectra are grouped to have a minimum of 40 counts per
# DESCRIPTION: group.
# DESCRIPTION: <p>
# DESCRIPTION: GIS RMFs are static and can are just retrieved from a collection
# DESCRIPTION: of pre-generated files.
# DESCRIPTION: SIS RMFs depend on epoch and location on the chip so they must be
# DESCRIPTION: generated using the FTOOL sisrmg.
# DESCRIPTION: The following proceedure is used to generate SIS RMFs in
# DESCRIPTION: 2 or 4 CCD mode:
# DESCRIPTION: <ol>
# DESCRIPTION: <li> Extract both a spectrum and an event list from the
# DESCRIPTION: extraction region.
# DESCRIPTION: <li> Divide the event list into an event list for each chip
# DESCRIPTION: and determine the number of counts on each chip.
# DESCRIPTION: on each chip.
# DESCRIPTION: <li> If the source is only on one chip, generate an RMF
# DESCRIPTION: using the spectrum generated above.
# DESCRIPTION: <li> If the source is spread over more than chip,
# DESCRIPTION: extract a spectrum from the event file for each chip
# DESCRIPTION: and then generate an RMF for each spectrum.
# DESCRIPTION: <li> sum the individual RMFs, weighting by the number of
# DESCRIPTION: events on each chip.
# DESCRIPTION: </ol>
# DESCRIPTION: <p>
# DESCRIPTION: An ARF is generated for each spectrum using
# DESCRIPTION: "extended source mode".
# DESCRIPTION: <p>
# DESCRIPTION: Blank sky spectra are extracted for each instrument and
# DESCRIPTION: source region file. The blank sky event files have been
# DESCRIPTION: previously selected to have COR (cut-off rigidity) consistent
# DESCRIPTION: with the selection criteria used to generate filtered
# DESCRIPTION: event files. The GIS blank sky event files have been
# DESCRIPTION: gisclean-ed. The spectra extracted from these blank sky
# DESCRIPTION: event files may be used for background subtraction
# DESCRIPTION: (with caution) in high galactic latitude sources.
# DESCRIPTION: <p>
# DESCRIPTION: Finally, each spectrum is plotted using xspec.
#
# VERSION: 2.1
#
# HISTORY: 0.0 -> 1.0 10/22/96
# HISTORY: fixed a bug where the source number was not incremented when
# HISTORY: the spectrum had fewer than the minimum number of counts
# HISTORY:
# HISTORY: 1.0 -> 1.1 12/19/96
# HISTORY: Modified in order to use the correct wmaps for sisrmg
# HISTORY: 1/3/97 Now handles MPC mode data
# HISTORY:
# HISTORY: 1.1 -> 1.2 1/22/96
# HISTORY: removed events${chip}.tmp in the case where there is more than one
# HISTORY: chip, but the counts are all on a single chip.
# HISTORY:
# HISTORY: 1.2 -> 1.3 2/6/97
# HISTORY: fixed the generate_filename call for rmf so that the source is
# HISTORY: in the source slot, not the index slot.
# HISTORY: 2/15/97 corrected problem with index variable and with RMF
# HISTORY: files not being cleare out before they are created.
# HISTORY: 2/17/97 corrected a bug where rmfs for individual chips were being
# HISTORY: generated using the spectra for the entire source instead of
# HISTORY: the individual chip.
# HISTORY: Also corrected a bug "$events${chip}.tmp" to "events${chip}.tmp"
# HISTORY: when removing event lists with no events. This was
# HISTORY: leaving behind temporary files.
# HISTORY: 2/18/97 Changed generate_filename call for GIS RMFs to accomodate
# HISTORY: different numbers of PHA bins.
# HISTORY: 2/19/97 xspec was hanging because of a mis-match in binning
# HISTORY: between SIS spectra and RMFs.
# HISTORY: This was caused by extracting spectra from
# HISTORY: individual chips with "group=no". Changed this to "group=yes".
# HISTORY: Filled out the description above.
# HISTORY:
# HISTORY: 1.3 -> 1.4 3/11/97
# HISTORY: put ignore bad into spectrum plotting to avoid excessive ink
# HISTORY: for leftover, ungrouped channels.
# HISTORY:
# HISTORY: 1.4 -> 1.5 3/25/97
# HISTORY: Added an explicit test for presence of sourcecat.
# HISTORY: 3/26/97
# HISTORY: Now used $STOOLS/group_event_files to group event files for
# HISTORY: spectra. As a result spectra with different active level
# HISTORY: discriminator settings will be kept separate.
# HISTORY: Also "commented out" the blanksky background spectra
# HISTORY: code by putting in an early "exit" until some better ideas
# HISTORY: come along for generating these.
# HISTORY:
# HISTORY: 1.5 -> 2.0
# HISTORY: Changed routine so that all extractions are done in
# HISTORY: detector coordinates
# HISTORY: Stopped generating source 0 spectra and responses unless there are
# HISTORY: no sources.
# HISTORY: Added an explicit check for whether a source was extracted.
# HISTORY: Added GIS deadtime correction
# HISTORY: Handles different GIS standard extraction regions for different
# HISTORY: spread discriminator settings.
# HISTORY:
# HISTORY: 2.0 -> 2.1 5/13/97
# HISTORY: Now increment the group index when there is no source catalog
#
##############################################################################
#DEBUG=1
$UTIL/milestone "Extracting spectra and generating response matricies"
############################
# some temporary file names
############################
xselcom=z.xco
listfile=listfile.tmp
#####################################################
# minimum counts per bin for grouping
#####################################################
minperbin=$($UTIL/read_parfile $PARFILE minperbin )
#########################################
# minimum events for keeping a PI file
#########################################
minpievents=$($UTIL/read_parfile $PARFILE minpievents )
#################################################################
# fetch some calibration files for response matrix generation
#################################################################
caldir=$( $UTIL/read_parfile $PARFILE caldir )
xrt_eff=$( $UTIL/read_parfile $PARFILE xrteff )
xrt_psf=$( $UTIL/read_parfile $PARFILE xrtpsf )
$UTIL/fetch_file $caldir $xrt_eff
$UTIL/fetch_file $caldir $xrt_psf
##############################
# region file source directory
##############################
regiondir=$($UTIL/read_parfile $PARFILE regiondir)
###########################
# SIS RMF calibration data
###########################
sispical=$($UTIL/read_parfile $PARFILE sispical )
echosh=$( $UTIL/read_parfile $PARFILE echosh )
rddhis=$( $UTIL/read_parfile $PARFILE rddhis )
$UTIL/fetch_file $caldir $sispical
$UTIL/fetch_file $caldir $echosh
$UTIL/fetch_file $caldir $rddhis
#############################################################################
# if there are multiple mkf files, then collect all the deatime info into
# one file
############################################################################
count=$($UTIL/read_parfile $JOBPAR frfcnt )
if [ $count -eq 1 ]; then
#########################################
# single filter file, just use that one
########################################
deadfile=$($UTIL/any_filename filter )
$UTIL/log_entry "Using $deadfile for GIS deadtime info"
else
###############################################
# combine multiple filter files
# note an alphabetical listing should give a
# the filter files in chronologivcal order.
###############################################
deadfile=deadfile.tmp
rm -f $deadfile
rm -f list.tmp
$UTIL/any_filename filter >list.tmp
$UTIL/log_entry "Collecting deadtime info from:"
$UTIL/file_to_log list.tmp
##########################################
# merge the files
##########################################
rm -f z.tmp
$UTIL/setup_parfile $FTOOLS/fmerge.par infiles="@list.tmp" \
outfile="z.tmp" \
columns="TIME G2_DEADT G3_DEADT"
$FTOOLS/fmerge >stdout.log 2>stderr.log
$UTIL/ftool_test fmerge $? $0 2 stdout.log stderr.log
rm -f list.tmp
################################
# sort them by time
################################
$UTIL/setup_parfile $FTOOLS/fmemsort.par infile="z.tmp" \
outfile=$deadfile \
column="TIME" \
method=insert \
ascend=yes \
copyall=no \
history=no
$FTOOLS/fmemsort >stdout.log 2>stdout.log
$UTIL/ftool_test fmemsort $? $0 2 stdout.log stderr.log
rm -f z.tmp
fi
#########################
# loop over instruments
#########################
for inst in s0 s1 g2 g3; do
index=101
case $inst in
s?)
modelist="02 12"
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
####################################################################
# create lists of event file to merge into spectra
#
####################################################################
#######################
# debugging output
########################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: mode=$mode
fi
###############################################
# group the event files
###############################################
rm -f allfiles.tmp
list=$($UTIL/generate_filename event $inst any $mode any )
ls $list >allfiles.tmp 2>/dev/null
if [ -s allfiles.tmp ]; then
##########################
# group the files
##########################
rm -f group.tmp
$STOOLS/group_event_files allfiles.tmp \
$LISTS/${full}_spectra_keys \
>stdout.log 2>stderr.log
code=$?
cp stdout.log group.tmp
$UTIL/stool_test group_event_files $code $0 2 \
stdout.log stderr.log
ngroups=$(tail -1 group.tmp |awk '{print $1}' )
else
#########################
# no files to group
#########################
ngroups=0
fi
rm -f allfiles.tmp
#######################
# debugging output
########################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: ngroups=$ngroups
echo ${0##*/}: list=$list
fi
##############################
# loop over event file groups
##############################
typeset -i group=1
while [ $group -le $ngroups ]; do
#########################################
# extract the event files in this group
#########################################
rm -f $listfile
awk '$1 == word {print $2}' word=$group group.tmp >$listfile
counts=$(awk '$1 == word {print $3}' word=$group group.tmp |\
head -1 )
if [ $counts -lt $minpievents ]; then
###############################################
# there aren't enough counts in the event file
###############################################
let group=$group+1
continue
fi
event=$(head -1 $listfile )
##############################################################
# determine image dimension to know which source list to use
##############################################################
case $inst in
s?) dimen=320;;
g?) dimen=$($UTIL/read_fits_keyword ${event}[0] RAWXBINS );;
*) $UTIL/exception $0 1 "Unknown instrument $inst";;
esac
####################################################
# determine a list of region files which mask off
# the "off-detector" regions
# These masks should not remove any events, but they
# are necessary to ensure that the WMAP is correct
####################################################
masklist=""
case "$inst" in
s?)
################################
# SIS,
# mask off the turned off chips
################################
ccdlist=$($UTIL/read_fits_keyword $event[0] \
S${inst#s}CCDLST)
chip=0
while [ chip -lt 4 ]; do
if [ -z "$(echo $ccdlist |grep $chip )" ]; then
########################
# this chip is off
########################
mask=$($UTIL/read_parfile $PARFILE \
${inst}notchip${chip} )
$UTIL/fetch_file $regiondir $mask
masklist="${masklist} $mask"
fi
let chip=$chip+1
done
#################################
# ... and the interchip gaps
#################################
mask=$($UTIL/read_parfile $PARFILE ${inst}offchip )
$UTIL/fetch_file $regiondir $mask
masklist="${masklist} $mask"
##################################
# and the area discriminated areas
##################################
teldef=$($UTIL/generate_filename teldef $inst )
discrim=$($UTIL/generate_filename discrim \
$inst $index $mode )
$STOOLS/discrim_mask $teldef $event >$discrim 2>stderr.log
$UTIL/stool_test discrim_mask $? $0 2 "" stderr.log
masklist="${masklist} $discrim"
;;
g?)
#########################################
# GIS
# mask off the outer ring and calsource
#########################################
spread=$($UTIL/read_fits_keyword ${event}[0] S_DSCR )
region=$($UTIL/read_parfile $PARFILE \
${inst}region${dimen}${spread})
$UTIL/fetch_file $regiondir $region
masklist="${masklist} $region"
;;
esac
#############################
# debugging output
#############################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: masklist=$masklist
fi
#########################
# loop over sources
#########################
sourcecat=$($UTIL/generate_filename sourcecat \
$full $dimen $sourcemode $bitrate )
###########################################################
# don't even make a background spectrum if there is no
# source catalog
##########################################################
if [ ! -a "$sourcecat" ]; then
$UTIL/log_entry "No source catalog $sourcecat"
let group=$group+1
continue
fi
if [ "$mode" = "71" -o "$mode" = "03" ]; then
####################
# MPC or FAST mode
####################
nsources=0
else
###############################################
# get the number of sources from the sourcecat
###############################################
nsources=$($UTIL/read_fits_keyword ${sourcecat}[1] NAXIS2 )
fi
#############################################################
# only create a source zero spectrum if there are no sources
#############################################################
if [ $nsources -eq 0 ]; then
typeset -i source=0
else
typeset -i source=1
fi
while [ $source -le $nsources ]; do
region=$($UTIL/generate_filename sourceregion $inst $dimen \
$sourcemode none \
$source)
spectrum=$($UTIL/generate_filename spectrum $inst \
$index $mode any $source)
xselcom=${spectrum%.*}.xco
rm -f $spectrum
$UTIL/log_entry "Extracting $spectrum from $region and:"
$UTIL/file_to_log $listfile
#######################
# debugging output
########################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: listfile=$listfile
echo ${0##*/}: index=$index
echo ${0##*/}: dimen=$dimen
echo ${0##*/}: sourcecat=$sourcecat
echo ${0##*/}: spectrum=$spectrum
fi
##################################
# 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 $(cat $listfile); do
echo "read events ${event}" >>$xselcom
done
echo "set phaname PI" >>$xselcom
if [ "$mode" != "71" -o "$mode" != "03" ]; then
#########################################
# no region filter for MPC or FAST modes
#########################################
for reg in $region $masklist; do
echo "filter region $reg" >>$xselcom
done
fi
echo "extract spectrum" >>$xselcom
echo "save spectrum ${spectrum} group=yes" >>$xselcom
case "$inst" in
s?) #################################################
# extract events as well for SIS multi-chip
# data. This is so the counts per/chip can be
# counted for proper weighting of the
# response matricies
################################################
ccdmod=$($UTIL/read_fits_keyword ${event}[0] \
S${inst#s}CCDMOD )
if [ $ccdmod -ne 1 ]; then
rm -f events.tmp
echo "extract events" >>$xselcom
echo "save events outfile=events.tmp use_events=yes"\
>>$xselcom
fi
;;
esac
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
#################################################
# check if spectrum was extracted OK
#################################################
if [ ! -a "$spectrum" ]; then
$UTIL/exception $0 2 "Spectrum $spectrum not extracted"
let source=$source+1
continue
fi
##########################################################
# if there were less than a minimum number of events,
# delete the PI file and skip to the next iteration
########################################################
totcts=$($UTIL/read_fits_keyword ${spectrum}[1] TOTCTS )
if [ $totcts -lt minpievents ]; then
$UTIL/log_entry \
"Deleting $spectrum since it has $totcts events"
rm -f $spectrum
rm -f events.tmp
let source=$source+1
continue
fi
##############################
# GIS deadtime correction
##############################
case "$inst" in
g?)
$UTIL/log_entry "Correcting $spectrum for dead time"
$UTIL/setup_parfile $FTOOLS/deadtime.par \
infile=$spectrum \
gtifile="$spectrum[2]" \
deadfile=$deadfile \
deadcol="G${inst#g}_DEADT" \
timecol=TIME \
gtistart=START \
gtistop=STOP \
exposure=ONTIME
$FTOOLS/deadtime >stdout.log 2>stderr.log
$UTIL/ftool_test deadtime $? $0 2 stdout.log stderr.log
;;
esac
#######################
# group the spectrum
#######################
$UTIL/log_entry \
"Grouping $spectrum to $minperbin counts/bin"
$UTIL/setup_parfile $FTOOLS/grppha.par \
infile=$spectrum \
outfile="!$spectrum" \
clobber=yes \
comm="group min $minperbin & show grouping & exit"
$FTOOLS/grppha >stdout.log 2> stderr.log
$UTIL/ftool_test grppha $? $0 2 stdout.log stderr.log
####################################
####################################
##
## generate response matricies
##
####################################
####################################
####################################################
# SIS and GIS instruments are treated differently
###################################################
case "$inst" in
s?) #######
# SIS
#######
#########################################
# how many CCDs does the source fall on?
#########################################
ccdmod=$($UTIL/read_fits_keyword ${spectrum}[0] \
S${inst#s}CCDMOD )
if [ $ccdmod -ne 1 ]; then
#############################
# multi-chip data
#############################
##########################################
# dump the chips that the extracted
# source events fall on
###########################################
rm -f z.tmp
$UTIL/dump_table "events.tmp[1]" CCDID >z.tmp
ccdlist=$(sort z.tmp |uniq)
ccdcount=$(echo $ccdlist |wc -w )
rm -f z.tmp
else
##############################################
# one CCD data
##############################################
ccdlist=$(echo $ccdlist |awk '{print $1}' )
ccdcount=1
fi
$UTIL/log_entry "Source occupies chip(s) $ccdlist"
###########################################
# generate RMF with a method depending on
# how many chips are spanned by the source
############################################
rmf=$($UTIL/generate_filename rmf $inst $index $mode \
any $source )
if [ $ccdcount -eq 1 ]; then
################################################
# source only spans one chip so just generate
# a single RMF
################################################
chip=$(echo $ccdlist)
###############################
# fetch correct ecddata file
##############################
split=$($UTIL/read_fits_keyword ${spectrum}[0] \
S${inst#s}_SPTR${chip} )
ecdata=$($UTIL/read_parfile $PARFILE \
ecd${inst#s}c${chip}p${split} )
$UTIL/fetch_file $caldir $ecdata
##############################
# generate RMF
##############################
$UTIL/log_entry "Generating $rmf"
rm -f $rmf
$UTIL/setup_parfile $FTOOLS/sisrmg.par \
infile=$spectrum \
arfile=NONE \
rmfile=$rmf \
datadir=./ \
ecdata=$ecdata \
phtopi=$sispical \
echosh=$echosh \
rddhis=$rddhis \
grades="0234" \
chip=$chip
$FTOOLS/sisrmg >stdout.log 2>stderr.log
$UTIL/ftool_test sisrmg $? $0 2 stdout.log \
stderr.log
else
####################################
# source spans multiple chips
####################################
###########################################
# extract a spectrum from the event file
# for each chip so that sisrmg will have
# a wmap for each chip
###########################################
for chip in $ccdlist; do
###########################################
# create a list of region files that mask
# off all but the current chip
###########################################
masklist=""
for chip2 in $ccdlist; do
if [ $chip -ne $chip2 ]; then
mask=$($UTIL/read_parfile $PARFILE\
${inst}notchip${chip} )
$UTIL/fetch_file $regiondir $mask
masklist="$masklist $mask"
fi
done
#########################
# extract the spectrum
#########################
$UTIL/log_entry \
"Extracting spectrum from chip $chip"
rm -f spec${chip}.tmp
rm -f $xselcom
echo "xsel" >$xselcom
echo "set mission ASCA" >>$xselcom
echo "set datadir ." >>$xselcom
echo "read events events.tmp" \
>>$xselcom
echo "set phaname PI" >>$xselcom
for reg in $masklist; do
echo "filter region $reg" >>$xselcom
done
echo "extract spectrum" >>$xselcom
echo "save spectrum spec${chip}.tmp group=yes" \
>>$xselcom
echo "exit save_session=no" >>$xselcom
$FTOOLS/xselect @${xselcom} </dev/null \
>stdout.log \
2>stderr.log
$UTIL/xselect_test $? $0 stdout.log \
stderr.log
rm -f $xselcom
#######################################
# generate an RMF for each active chip
#######################################
###############################
# fetch correct ecddata file
##############################
split=$($UTIL/read_fits_keyword \
${spectrum}[0] \
S${inst#s}_SPTR${chip} )
ecdata=$($UTIL/read_parfile $PARFILE \
ecd${inst#s}c${chip}p${split} )
$UTIL/fetch_file $caldir $ecdata
########################
# generate RMF
########################
$UTIL/log_entry \
"Generating RMF for chip $chip"
rm -f rmf${chip}.tmp
$UTIL/setup_parfile $FTOOLS/sisrmg.par \
infile=spec${chip}.tmp \
arfile=NONE \
rmfile=rmf${chip}.tmp \
datadir=./ \
ecdata=$ecdata \
phtopi=$sispical \
echosh=$echosh \
rddhis=$rddhis \
grades="0234" \
chip=$chip
$FTOOLS/sisrmg >stdout.log 2>stderr.log
$UTIL/ftool_test sisrmg $? $0 2 stdout.log \
stderr.log
done
#######################################
# add the RMFs weighted by counts/chip
#######################################
rm -f wgtlist.tmp
for chip in $ccdlist; do
nevents=$($UTIL/read_fits_keyword \
spec${chip}.tmp[0] \
TOTCTS )
wgt=$($STOOLS/equals \
"${nevents}/${totcts}" )
echo rmf${chip}.tmp $wgt >>wgtlist.tmp
rm -f spec${chip}.tmp
done
$UTIL/log_entry "Chip weights:"
$UTIL/file_to_log wgtlist.tmp
rm -f $rmf
$UTIL/setup_parfile $FTOOLS/addrmf.par \
listfile=wgtlist.tmp \
rmffile=$rmf \
clobber=yes
$FTOOLS/addrmf >stdout.log 2>stderr.log
$UTIL/ftool_test addrmf $? $0 2 stdout.log \
stderr.log
rm -f rmf?.tmp
rm -f wgtlist.tmp
fi
################################################
# dummy values for GIS ARF calibration files
###############################################
bethick=NONE
grid=NONE
;;
g?) ######################################
# GIS -- RMF is fetched, not generated
######################################
rmf=$($UTIL/generate_filename rmf \
$inst $index $mode any $source)
$UTIL/fetch_file $caldir $rmf
############################
# GIS ARF calibration files
############################
bethick=$($UTIL/read_parfile $PARFILE ${inst}bethick )
grid=$( $UTIL/read_parfile $PARFILE ${inst}grid )
$UTIL/fetch_file $caldir $bethick
$UTIL/fetch_file $caldir $grid
;;
*) $UTIL/exception $0 1 "Unknown instrument $inst";;
esac
##########################
# generate an ARF
##########################
arf=$($UTIL/generate_filename arf $inst \
$index $mode any $source)
############################
# extended or compact?
############################
if [ $source -eq 0 ]; then
#######################################
# background region is always extended
#######################################
point=no
else
##############################
# check the sourcecat
##############################
ext=$($UTIL/read_bintable $sourcecat extended $source)
case "$ext" in
T) point=no;;
F) point=yes;;
*) $UTIL/exception $0 2 "Unknown extended flag: $ext";;
esac
fi
$UTIL/log_entry "Generating $arf with point=$point"
$UTIL/setup_parfile $FTOOLS/ascaarf.par \
phafile=${spectrum}\
rmffile=${rmf} \
outfile=${arf} \
point=$point \
simple=yes \
xrtrsp=${xrt_eff} \
xrtpsf=${xrt_psf} \
bethick=${bethick} \
grid=${grid} \
$FTOOLS/ascaarf >stdout.log 2>stderr.log
$UTIL/ftool_test ascaarf $? $0 2 stdout.log stderr.log
#################################################
# tell the spectrum file the names of the
# response matrix files
####################################################
$UTIL/add_fits_keyword ${spectrum}[1] RESPFILE "'${rmf}'"
$UTIL/add_fits_keyword ${spectrum}[1] ANCRFILE "'${arf}'"
let source=$source+1
done #loop over sources
rm -f $listfile
let index=$index+1
let group=$group+1
done #loop over event file groups
rm -f group.tmp
done #loop over mode
done #loop over instrument
############################################
# remove merged deadtime file
############################################
count=$($UTIL/read_parfile $JOBPAR frfcnt )
if [ $count -gt 1 ]; then
##########################
# there are multiple frfs
##########################
rm -f $deadfile
fi
#####################################################################
#####################################################################
##
## background spectra from blanksky observations
##
#####################################################################
######################################################################
for inst in s0 s1 g2 g3; do
###################################################################
# don't extract background spectra - at least for the time being
##################################################################
continue
$UTIL/log_entry "Creating $inst blanksky backgrounds"
case "$inst" in
s*) mode=02
full=sis
dimen=320;;
g*) mode=70
full=gis
dimen=256;;
*) $UTIL/exception $0 1 "Unknown instrument $inst";;
esac
sourcecat=$($UTIL/generate_filename sourcecat $full $dimen $mode any)
if [ ! -a "$sourcecat" ]; then
$UTIL/log_entry "$sourcecat does not exist"
continue
fi
################################
# get the blank sky event file
###############################
blanksky=$($UTIL/read_parfile $PARFILE ${inst}blanksky )
caldir=$( $UTIL/read_parfile $PARFILE caldir )
$UTIL/fetch_file $caldir $blanksky
#########################
# loop over sources
#########################
nsources=$($UTIL/read_fits_keyword ${sourcecat}[1] NAXIS2 )
typeset -i source=0
while [ $source -le $nsources ]; do
backblank=$($UTIL/generate_filename backblank $inst $dimen \
$mode any $source)
rm -f $backblank
$UTIL/log_entry "Extracting $backblank"
region=$($UTIL/generate_filename sourceregion $inst $dimen \
$mode none $source)
##################################
# generate an xselect command file
##################################
rm -f $xselcom
echo "xsel" >$xselcom
# echo "%ED%off" >>$xselcom
echo "set mission ASCA" >>$xselcom
echo "set datadir ." >>$xselcom
echo "read events $blanksky" >>$xselcom
echo "set phaname PI" >>$xselcom
echo "filter region \"$region\"" >>$xselcom
echo "extract spectrum" >>$xselcom
echo "save spectrum ${backblank} group=yes" >>$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
let source=$source+1
done
done
###############################################
# plot all of the spectra
##############################################
XANADU=$($UTIL/read_parfile $PARFILE xanadu )
list=$($UTIL/generate_filename spectrum any any any any any ) list=$(ls $list 2>/dev/null)
command=z.xcm
for spectrum in $list; do
inst=$( $UTIL/parse_filename inst $spectrum )
index=$( $UTIL/parse_filename index $spectrum )
mode=$( $UTIL/parse_filename mode $spectrum )
source=$($UTIL/parse_filename source $spectrum )
specplot=$($UTIL/generate_filename specplot $inst $index $mode \
any $source )
#pibins=$($UTIL/read_fits_keyword $spectrum[1] NAXIS2 )
#rmf=$($UTIL/read_fits_keyword $spectrum[1] RESPFILE )
#if [ ! -r "$rmf" ]; then
# $UTIL/exception $0 2 \
# "Can't read RMF $rmf for $spectrum, skipping plot"
# continue
#fi
#rmfbins=$($UTIL/read_fits_keyword $rmf[1] NAXIS2 )
#
#if [ "$pibins" != "$rmfbins" ]; then
# $UTIL/exception $0 2 \
# "$pibins in $spectrum and $rmfbins in $rmf, skipping plot"
# continue
#fi
$UTIL/log_entry "Plotting $specplot from $spectrum"
rm -f $command
echo "data $spectrum" >>$command
echo "setplot energy" >>$command
echo "ignore bad" >>$command
echo "iplot" >>$command
echo "log y on" >>$command
echo "rescale y" >>$command
echo "hard /ps" >>$command
echo "quit" >>$command
echo "exit" >>$command
echo "y" >>$command
$XANADU/sol/bin/xspec @$command </dev/null >stdout.log 2>stderr.log
$UTIL/stool_test xspec $? $0 2 stdout.log stderr.log
mv -f pgplot.ps $specplot
rm -f $command
done
exit 0