calsource (version 1.0)
You can also look at:
###########################################################################
#
#
# SYNTAX: calsource
#
# BRIEFLY: Extract GIS calibration source spectra
#
# DESCRIPTION: This routine extracts spectra of the GIS calibration source
#
# VERSION: 1.0
#
# HISTORY: 0.0 -> 1.0 5/12/97
# HISTPRY: added missing explicit path to exception
#
##############################################################################
#DEBUG=1
command=z.xco
regiondir=$($UTIL/read_parfile $PARFILE regiondir )
caldir=$( $UTIL/read_parfile $PARFILE caldir )
#####################################################
# minimum counts per bin for grouping
#####################################################
#minperbin=$($UTIL/read_parfile $PARFILE minperbin )
calpeaktoll=$($UTIL/read_parfile $PARFILE calpeaktoll)
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
$STOOLS/group_event_files files.tmp $LISTS/keywords_calsource |\
grep "^[0-9]" >sorted.tmp
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 "filter region $region" >>$command
echo "extract spectrum" >>$command
echo "save spectrum ${calsource} group=yes" >>$command
echo "exit save_session=no" >>$command
#################
# run xselect
#################
$FTOOLS/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
continue
fi
#############################
# 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 $FTOOLS/grppha.par \
# infile=$calsource \
# outfile="!$calsource" \
# 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
##############################
# 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 **-5.0" >>$command
echo "ignore 7.0-**" >>$command
echo "model const/b + gauss/b + gauss/b " >>$command
echo " " >>$command # constant 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 5 = 2 1.1" >>$command # peak energies
#echo "newpar 7 = 4 .2 " >>$command # line ratio
echo "newpar 3 = 6" >>$command # both lines have the same width
echo "renorm" >>$command
echo "fit" >>$command
#echo "error 2" >>$command
echo "setplot energy" >>$command
echo "iplot data" >>$command
echo "label top GIS${inst#g} Calibration Source," \
"should peak near 5.9 keV" >>$command
echo "log x off" >>$command
echo "rescale x 5 7" >>$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
code=$?
#######################################
# parse the sdtout to get fit results
#######################################
best=$(awk '/XSPEC> fit/, /XSPEC> error/' stdout.log |\
grep LineE |head -1 |awk '{print $6}' )
$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 spectrum $calsource peaks at $best keV"
error=$($STOOLS/equals "abs($best - 5.9 )" )
thing=$($STOOLS/floatcmp $error $calpeaktoll )
if [ "$thing" = "gt" ]; then
$UTIL/exception $0 2 "GIS${inst#g} gain problem in $calsource"
fi
################################################################
# 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)"
fi
let group=$group+1
done
rm -f sorted.tmp
done
exit 0