calsource (version 2.4)
You can also look at:
#! /usr/bin/ksh
###########################################################################
#
# SYNTAX: calsource
#
# BRIEFLY: Extract GIS calibration source spectra
#
# DESCRIPTION: This routine extracts spectra of the GIS calibration sources.
# DESCRIPTION: It then uses <TT>xspec</TT> to plot these spectra and find their
# DESCRIPTION: peak. The calibration source spectrum is fitted with
# DESCRIPTION: a pair of Gaussians whose peak energies are constrained
# DESCRIPTION: to have the ratio of 1.1 and whose widths are constrained
# DESCRIPTION: to be identical. The two peaks are needed to represent the
# DESCRIPTION: K alpha and K beta peaks from the calibration source.
# DESCRIPTION: If the lower energy (K alpha) peak does not fall within
# DESCRIPTION: %calpeaktoll keV of 5.9 keV, then an error is given.
#
# VERSION: 2.4
#
# HISTORY: 0.0 -> 1.0 5/12/97
# HISTORY: Added missing explicit path to exception.
# HISTORY:
# HISTORY: 1.0 -> 1.1 5/23/97
# HISTORY: Now runs stool test after group_event_files run.
# HISTORY: Fixed a bug where the index variable was not incremented.
# HISTORY: Now, if there are multiple calsource spectra for a given instrument,
# HISTORY: only issues an error if the "bestgroup" (i.e. the one with the
# HISTORY: most events in the unfiltered files) has a cal. peak outside
# HISTORY: of the allowed tolerance.
# HISTORY:
# HISTORY: 1.1 -> 1.2 8/4/97
# HISTORY: Now skips if there are no .unf files for a given instrument.
# HISTORY:
# HISTORY: 1.2 -> 2.0 8/7/97
# HISTORY: Made some slight changes to the fitting proceedure to match what
# HISTORY: Ken Ebisawa reports the GIS team uses. In one test case the changes
# HISTORY: accounted for a difference of only few x .001 keV in the peak
# HISTORY: energy.
# HISTORY: Also started extracting the fit error from the log and storing
# HISTORY: it in the job.par.
# HISTORY:
# HISTORY: 2.0 -> 2.1 1998-09-17
# HISTORY: modified to use xanadu initialization
# HISTORY:
# HISTORY: 2.1 -> 2.2 1999-07-07
# HISTORY: Supplied missing script designation at the top - not sure why this
# HISTORY: worked before.
# HISTORY:
# HISTORY: 2.2 -> 2.3 1999-11-30
# HISTORY: changed to use calpeaklo and calpeakhi instead of the single
# HISTORY: calpeaktoll value.
# HISTORY:
# HISTORY: 2.3 -> 2.4 2000-02-02
# HISTORY: Updated xspec syntax from "xspec @command" to "xspec - command"
# HISTORY: and modified command script to remove "iplot".
# HISTORY: Also added code to work around an extractor bug in WMAPS.
# HISTORY: Modified xspec fitting script to use a powerlaw with index 0.0
# HISTORY: instead of "const" to model the continuum. Const is a
# HISTORY: multiplicative model, though older versions of xspec
# HISTORY: apparently allowed it ot be used as an additive model.
# HISTORY: Changed xspec log parsing to accomodate new xspec version.
# HISTORY: Now explicitly specify detector coordinates for xselect
#
###########################################################################
#DEBUG=1
command=z.xco
regiondir=$($UTIL/read_parfile $PARFILE regiondir)
caldir=$( $UTIL/read_parfile $PARFILE caldir)
######################################
# Minimum counts per bin for grouping
######################################
calpeaklo=$($UTIL/read_parfile $PARFILE calpeaklo)
calpeakhi=$($UTIL/read_parfile $PARFILE calpeakhi)
mode=70;
$UTIL/milestone "Extracting GIS calibration source spectra"
for inst in g2 g3; do
index=201;
rm -f files.tmp
rm -f sorted.tmp
list=$($UTIL/generate_filename unfiltered $inst any $mode any)
ls $list >files.tmp 2>/dev/null
if [ ! -s files.tmp ]; then
###########################################
# No files, so skip to the next instrument
###########################################
$UTIL/log_entry "No unfiltered $inst files"
continue
fi
$STOOL/group_event_files files.tmp $LISTS/keywords_calsource \
>stdout.log 2>stderr.log
code=$?
cp stdout.log sorted.tmp
$UTIL/stool_test group_event_files $code $0 2 stdout.log stderr.log
rm -f files.tmp
ngroups=$(tail -1 sorted.tmp | awk '{print $1}')
bestgroup=$(awk '{print $3" "$1}' sorted.tmp | uniq | sort | tail -1 | \
awk '{print $2}')
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: ngroups=$ngroups
echo ${0##*/}: bestgroup=$bestgroup
fi
typeset -i group=1
while [ $group -le $ngroups ]; do
list=$( awk '$1 == word {print $2}' word=$group sorted.tmp)
first=$( awk '$1 == word {print $2}' word=$group sorted.tmp | head -1)
nevents=$(awk '$1 == word {print $3}' word=$group sorted.tmp | head -1)
rawxbins=$($UTIL/read_fits_keyword $first[0] RAWXBINS)
# phabins=$($UTIL/read_fits_keyword $first[0] PHABINS)
##########################
# Fetch the region filter
##########################
region=$($UTIL/read_parfile $PARFILE ${inst}calregion${rawxbins})
$UTIL/fetch_file $regiondir $region
#####################
# Spectrum file name
#####################
calsource=$($UTIL/generate_filename calsource $inst $index $mode)
rm -f $calsource
$UTIL/log_entry "Extracting $calsource from $list"
#######################
# xselect command file
#######################
rm -f $command
echo "xsel" >>$command
echo "set mission ASCA" >>$command
echo "set datadir ." >>$command
for file in $list; do
echo "read events $file " >>$command
done
echo "set image det" >>$command
echo "filter region $region" >>$command
echo "extract spectrum" >>$command
echo "save spectrum ${calsource} group=yes" >>$command
echo "exit save_session=no" >>$command
##############
# Run xselect
##############
$FTOOL/xselect @${command} </dev/null > stdout.log 2> stderr.log
$UTIL/xselect_test $? $0 stdout.log stderr.log
rm -f $command
rm -f list1.tmp
rm -f list2.tmp
############################################
# Skip the rest if no spectrum was produced
############################################
if [ ! -s "$calsource" ]; then
rm -f $calsource
let group=$group+1
continue
fi
##################################
# fix an extractor bug in spectra
##################################
$UTIL/unblank_wmap $calsource
##############################
# Add name of RMF to spectrum
##############################
rmf=$($UTIL/generate_filename rmf $inst $index $mode)
$UTIL/fetch_file $caldir $rmf
$UTIL/add_fits_keyword ${calsource}[1] RESPFILE "'${rmf}'"
#####################
# Group the spectrum
#####################
#$UTIL/log_entry "Grouping $calsource to $minperbin counts/bin"
#$UTIL/setup_parfile $FTOOL/grppha.par \
# infile=$calsource \
# outfile="!$calsource" \
# 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
##############################
# Plot spectrum and fit model
##############################
plot=$($UTIL/generate_filename calsourceplot $inst $index $mode)
$UTIL/log_entry "Plotting $plot and finding peak"
rm -f $command
echo "query yes" >>$command
echo "data $calsource" >>$command
echo "ignore bad" >>$command
echo "ignore **-4.0" >>$command
echo "ignore 8.0-**" >>$command
echo "model powerlaw/b + gauss/b + gauss/b " >>$command
echo "0.0" >>$command # continuum slope
echo " " >>$command # continuum normalization
echo "5.9" >>$command # initial K alpha line energy in keV
echo "0.0" >>$command # delta function
echo " " >>$command # normalization
echo "6.5" >>$command # initial K beta line energy in keV
echo "0.0" >>$command # delta function
echo " " >>$command # normalization
echo "newpar 6 = 3 1.101005" >>$command # peak energies
echo "newpar 7 = 4 1.049288" >>$command #line widths scale ~ sqrt(E)
echo "freeze 1" >>$command # continuum slope
echo "renorm" >>$command
echo "fit" >>$command
#echo "error 2" >>$command
echo "setplot energy" >>$command
echo "setplot com label top GIS${inst#g} Calibration Source," \
"should peak near 5.9 keV" >>$command
echo "setplot com log x off" >>$command
echo "setplot com rescale x 4 8" >>$command
echo "setplot com hard /ps" >>$command
echo "plot" >>$command
echo "exit" >>$command
$XANBIN/bin/xspec - $command </dev/null >stdout.log 2>stderr.log
code=$?
rm -f xautosav.xcm
######################################
# Parse the sdtout to get fit results
######################################
awk '/XSPEC> fit/, /XSPEC> error/' stdout.log | \
grep LineE | head -1 | \
read dum1 dum2 dum3 dum4 dum5 dum6 best dum7 errorbar
errorbar=$($STOOL/equals $errorbar)
$UTIL/stool_test xspec $code $0 2 stdout.log stderr.log
mv -f pgplot.ps $plot
rm -f $command
####################################################
# Repeat the peak energy and check if it is correct
####################################################
$UTIL/log_entry \
"Calibration source $calsource peaks at $best +/- $errorbar keV"
###############################################################
# If this is the best group, store the calpeak in the par file
###############################################################
if [ $group -eq $bestgroup ]; then
$UTIL/add_parameter $JOBPAR ${inst}calpeak "$best" r \
"GIS${inst#g} calibration source peak (keV)"
$UTIL/add_parameter $JOBPAR ${inst}calerr "$errorbar" r \
"GIS${inst#g} calibration source fit error (keV)"
thinglo=$($STOOL/floatcmp $best $calpeaklo )
thinghi=$($STOOL/floatcmp $best $calpeakhi )
if [ "$thinglog" = "lt" -o $thinghi = "gt" ]; then
$UTIL/exception $0 2 \
"GIS${inst#g} gain problem in $calsource"
fi
fi
let index=$index+1
let group=$group+1
done
rm -f sorted.tmp
done
exit 0