spectra (version 3.2)
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: matrices, and 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
# DESCRIPTION: %minperbin counts per group, and
# DESCRIPTION: spectra are not extracted if they would have fewer than
# DESCRIPTION: %minpievents events.
# DESCRIPTION: GIS spectra are corrected for dead time using the
# DESCRIPTION: <TT>deadtime</TT> FTOOL.
# DESCRIPTION: <P>
# DESCRIPTION: FITS spectrum files contain a WMAP extension containing
# DESCRIPTION: an image of the extraction region with pixels outside the
# DESCRIPTION: region set to -1. The FTOOL <TT>ascaarf</TT> uses this
# DESCRIPTION: extension to determine what parts of the focal plane
# DESCRIPTION: should be considered when calculating the ancillary
# DESCRIPTION: response matrix. Since event files do not carry any information
# DESCRIPTION: about previous region filters or the field of view of the
# DESCRIPTION: detectors, these things must be explicitly accounted for when
# DESCRIPTION: extracting spectra. The source region filter is followed
# DESCRIPTION: by a number of other region filters which mask off
# DESCRIPTION: "forbidden" regions of the focal plane. For the GIS this
# DESCRIPTION: means using the standard region filter, which excludes
# DESCRIPTION: the noisy outer ring and calibration source. For the SIS, things
# DESCRIPTION: are a little more complicated. A set of region filters are
# DESCRIPTION: used to mask off the inter-chip gaps and inactive CCD chips.
# DESCRIPTION: Also, a set of region filters are generated to mask off the
# DESCRIPTION: parts of the chips which are excluded by the area discriminators
# DESCRIPTION: onboard ASCA.
# DESCRIPTION: <P>
# DESCRIPTION: GIS RMFs are static, and are retrieved from a collection
# DESCRIPTION: of pre-generated files.
# DESCRIPTION: SIS RMFs depend on epoch and location on the chip,
# DESCRIPTION: so they must be
# DESCRIPTION: generated using the FTOOL <TT>sisrmg</TT>.
# DESCRIPTION: The following procedure 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: <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: the <TT>ascaarf</TT> FTOOL.
# DESCRIPTION: <P>
# DESCRIPTION: Finally, each spectrum is plotted using xspec.
#
# VERSION: 3.2
#
# 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 cleared 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 uses $STOOL/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 increments the group index when there is no source catalog.
# HISTORY:
# HISTORY: 2.1 -> 2.2 5/14/97
# HISTORY: Added more explicit checking for valid RMFs and ARFs.
# HISTORY: Accomodated new parameter syntax in FTOOLS 4.0 addrmf,
# HISTORY: though it is still backwards compatible with FTOOLS 3.x
# HISTORY: Corrected a bug which used "chip" instead of "chip2" when selecting
# HISTORY: masks to extract spectra from individual chips to make RMFs.
# HISTORY:
# HISTORY: 2.2 -> 2.3 6/0/97
# HISTORY: Added some "rm -f $listfile"s.
# HISTORY:
# HISTORY: 2.3 -> 2.4 7/9/97
# HISTORY: Now skips spectra for BRIGHT mode files with split threshold
# HISTORY: not equal to 40. This is becuase the ecd calibration files
# HISTORY: for SPLIT=20 do not cover bright mode. The theory (according to
# HISTORY: Geoff Crew) is that people would only want split=20 to get higher
# HISTORY: spectral resolution, in which case they will also use BRIGHT2 mode.
# HISTORY:
# HISTORY: 2.4 -> 2.5 7/25/97
# HISTORY: Fixed a bug in the loop over chips which caused SIS region
# HISTORY: filter masks to be chosen improperly. Also, the masklist variable
# HISTORY: was being reused when selecting individual chips to generate RMFS.
# HISTORY:
# HISTORY: 2.5 -> 2.6 8/6/97
# HISTORY: Added a check for whether or not a source region file exists
# HISTORY: before trying to extract with it.
# HISTORY:
# HISTORY: 2.6 -> 2.7 9/25/97
# HISTORY: Before FTOOLS 4.0, the deadtime FTOOL was not able to handle
# HISTORY: multiple mkf files, so this routine had code to merge all
# HISTORY: mkf files. This code was removed and a list of all mkf files
# HISTORY: is now given to deadtime.
# HISTORY:
# HISTORY: 2.7 -> 2.8 10/14/97
# HISTORY: Deleted the code which removed the merged deadtime file. This
# HISTORY: was obsolete since the last change.
# HISTORY:
# HISTORY: 2.8 -> 2.9 4/16/98
# HISTORY: No longer give the "spectrum not extracted" warning since
# HISTORY: this is more of an annoyance then a help.
# HISTORY:
# HISTORY: 2.9 -> 3.0 1998-09-17
# HISTORY: modified to use cxanadu initialization
# HISTORY:
# HISTORY: 3.0 -> 3.1 1999-03-24
# HISTORY: Deleted the check for pre-FTOOLS 4.0 syntax for addrmf.
# HISTORY: Streamlined code which determines "ccdlist".
# HISTORY:
# HISTORY: 3.1 -> 3.2 2000-02-02
# HISTORY: Changed xspec syntax from "xspec @command" to "xspec - command"
# HISTORY: and modified command syntax since iplot is no longer allowed.
# HISTORY: Also added code to fix an extractor bug in WMAPS.
# HISTORY: Now explicitly specify detector coordinates for xselect.
# HISTORY: changed the exposure parameter in deadtime from "ONTIME" to
# HISTORY: "EXPOSURE".
#
############################################################################
#DEBUG=1
$UTIL/milestone "Extracting spectra and generating response matrices"
############################
# 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
########################
# 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 files 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
$STOOL/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
rm -f $listfile
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
########################
# Check split threshold
########################
if [ "$mode" = "02" ]; then
#################################################
# Do not extract BRIGHT mode spectra if the
# split threshold = 20
# This is because the ecd calibration file
# needed to make an RMF doesn't cover this case.
#################################################
split=$($UTIL/read_fits_keyword $event[0] S${inst#s}_SPTR0)
if [ "$split" != "40" ]; then
$UTIL/log_entry \
"Skipping BRIGHT mode files with SPTR = $split"
$UTIL/file_to_log $listfile
rm -f $listfile
let group=$group+1
continue
fi
fi
#####################################################
# 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)
for chip in 0 1 2 3; 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
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)
$STOOL/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
rm -f $listfile
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)
##############################
# Check if region file exists
##############################
if [ ! -a "$region" ]; then
$UTIL/log_entry \
"Skipping $spectrum since $region does not exist"
let source=$source+1
continue
fi
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
echo "set image det" >>$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 matrices
##############################################
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 command file to xselect
###################################
$FTOOL/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
######################################
# fix an extractor bug in FTOOLS 5.0
######################################
$UTIL/unblank_wmap $spectrum
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: totcts=$totcts
fi
##########################
# GIS deadtime correction
##########################
case "$inst" in
g?)
$UTIL/log_entry "Correcting $spectrum for dead time"
rm -f deadlist.tmp
ls ft*.mkf >deadlist.tmp
$UTIL/setup_parfile $FTOOL/deadtime.par \
infile=$spectrum \
gtifile="$spectrum[2]" \
deadfile=@deadlist.tmp \
deadcol="G${inst#g}_DEADT" \
timecol=TIME \
gtistart=START \
gtistop=STOP \
exposure=EXPOSURE
$FTOOL/deadtime >stdout.log 2>stderr.log
$UTIL/ftool_test deadtime $? $0 2 stdout.log stderr.log
rm -f deadlist.tmp
;;
esac
#####################
# Group the spectrum
#####################
$UTIL/log_entry \
"Grouping $spectrum to $minperbin counts/bin"
$UTIL/setup_parfile $FTOOL/grppha.par \
infile=$spectrum \
outfile="!$spectrum" \
clobber=yes \
comm="group min $minperbin & show grouping & exit"
$FTOOL/grppha >stdout.log 2> stderr.log
$UTIL/ftool_test grppha $? $0 2 stdout.log stderr.log
##############################
##############################
##
## Generate response matrices
##
##############################
##############################
##################################################
# 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
####################################
ccdlist=$($UTIL/dump_table "events.tmp[1]" CCDID | \
sort -u)
ccdcount=$(echo $ccdlist | wc -w)
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 $FTOOL/sisrmg.par \
infile=$spectrum \
arfile=NONE \
rmfile=$rmf \
datadir=./ \
ecdata=$ecdata \
phtopi=$sispical \
echosh=$echosh \
rddhis=$rddhis \
grades="0234" \
chip=$chip
$FTOOL/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
# note we only mask off chips from which
# there are counts. This carries with it the
# assumption that if the source extraction
# region falls on a given chip then there
# will be at least one count from that chip
############################################
chipmasklist=""
for chip2 in $ccdlist; do
if [ $chip -ne $chip2 ]; then
mask=$($UTIL/read_parfile $PARFILE \
${inst}notchip${chip2})
$UTIL/fetch_file $regiondir $mask
chipmasklist="$chipmasklist $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
echo "set image det" >>$xselcom
for reg in $region \
$masklist \
$chipmasklist; 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
$FTOOL/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 $FTOOL/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
$FTOOL/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=$($STOOL/equals "${nevents}/${totcts}")
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: chip=$chip
echo ${0##*/}: nevents=$nevents
echo ${0##*/}: totcts=$totcts
fi
echo rmf${chip}.tmp $wgt >>wgtlist.tmp
rm -f spec${chip}.tmp
done
$UTIL/log_entry "Chip weights:"
$UTIL/file_to_log wgtlist.tmp
###############
# Add the rmfs
###############
rm -f $rmf
$UTIL/setup_parfile $FTOOL/addrmf.par \
list=@wgtlist.tmp \
rmffile=$rmf \
clobber=yes
$FTOOL/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 $FTOOL/ascaarf.par \
phafile=${spectrum} \
rmffile=${rmf} \
outfile=${arf} \
point=$point \
simple=yes \
xrtrsp=${xrt_eff} \
xrtpsf=${xrt_psf} \
bethick=${bethick} \
grid=${grid}
$FTOOL/ascaarf >stdout.log 2>stderr.log
$UTIL/ftool_test ascaarf $? $0 2 stdout.log stderr.log
if [ $? -ne 0 ]; then
exit 1
fi
##########################################
# Tell the spectrum file the names of the
# response matrix files
##########################################
if [ -s "$rmf" ]; then
$UTIL/add_fits_keyword ${spectrum}[1] RESPFILE "'${rmf}'"
else
$UTIL/exception $0 2 "Couldn't generate $rmf"
rm -f $rmf
fi
if [ -s "$arf" ]; then
$UTIL/add_fits_keyword ${spectrum}[1] ANCRFILE "'${arf}'"
else
$UTIL/exception $0 2 "Couldn't generate $arf"
rm -f $arf
fi
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
##################################################
##################################################
##
## 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 "set image det" >>$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
##################################
$FTOOL/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
##########################
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)
##########################
# Check response matrices
##########################
rmf=$($UTIL/read_fits_keyword ${spectrum}[1] RESPFILE)
arf=$($UTIL/read_fits_keyword ${spectrum}[1] ANCRFILE)
if [ ! -s "$rmf" -o ! -s "$arf" ]; then
###########################################
# Can't read response matrix, so skip plot
###########################################
$UTIL/log_entry "Can't read $arf or $rmf, skipping $specplot"
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 "setplot com log y on" >>$command
echo "setplot com rescale y" >>$command
echo "setplot com hard /ps" >>$command
echo "plot" >>$command
echo "exit" >>$command
# echo "y" >>$command
$XANBIN/bin/xspec - $command </dev/null >stdout.log 2>stderr.log
$UTIL/stool_test xspec $? $0 2 stdout.log stderr.log
rm -f xautosav.xcm
mv -f pgplot.ps $specplot
rm -f $command
done
exit 0