Initial script

#!/star/local/bin/tcsh 
nice +19
######################################################################
#  initial.csh RAWDIR ROR ROOTNAME BGRADIN BGRADOUT RA DEC
#
#  preliminary source search and image creation. Background
#  subtraction will be rerun until no more new sources are found.
#  If back ground needs scaling this should be done automatically.
#  Will then produce background subtracted corrected spectral cube, 
#  background model spectral cube and some images centered on the 
#  user supplied group centre.
#
#  Based loosely on dks version of pdjb's image1.csh and image2.csh
#  Stephen Helsdon
#  Last modified 9/12/99  (VERSION 3.0b3)
#
#  1.2 - added radius variable to xrtcorr (now removed)
#  1.3 - corrected error in call to getspec (now removed)
#  1.4 - added saving of xpos and ypos. Added calls to extra.csh
#        remmoved rubbish comments
#  2.0 - restructured slightly, changed to new version of wellard
#  2.1 - temporary version - working on setting up for new reductions
#        will eventually result in V3 when finished.
#  3.0b1 - beta version of script - major changes - changed command 
#        syntax to fit current version of asterix (2.3b1), improved
#        several parts of the script (does things more neatly), added
#        extra stuff such as pkstart to start up necessary packages, 
#        default.set to set up some basic defaults which can be 
#        changed (adds a bit of versatility, and shouldn't now cause
#        the rest of the script to crash as happened previously if 
#        you changed values within it). The script also now tries to 
#        be aware of what type of computer you are running on as this 
#        has some effect on exactly what is run. Also a little more 
#        ready (hopefully) to adapt to different conditions. Seems to 
#        work ok on suns (early tests anyway). Should work ok on alphas
#        but not yet tested. Not yet working on linux as stacks of the
#        installed software is not working properly...
#  3.0b2 - changed format of extra call...
#  3.0b3 - further small changes - needs people to test it now!
#  3.0b4 - couple of small changes - Linux now same as for other 
#        machines.
#####################################################################
#
if ( $#argv != 7 ) then
  echo "Usage: initial.csh  dir  ror  rootname  inrad outrad  RA  DEC"
  echo "See initial.txt for more details"
  exit
endif            
# start asterix, obelix and kappa
source pkstart
#
#####################################################################
#
echo 'INITIAL.CSH VERSION 3.0b4'
#  create an empty source list file (needed later)
touch old_ard.dat
#
set itt=1
#
# set up global variables
# 
set rawdir=$1
set ror=$2
set root=$3
set bgradin=$4
set bgradout=$5
set ra=$6
set dec=$7
#
# set up some defaults if none exist or read in current values from 
# defaults file
source default.set
#
#####################################################################
#
#  uncompress data (only one of uncompress or gunzip will be needed)
gunzip -dq *.gz
if (-e $ror) tar xvf $ror
if (-e $ror'.tar') tar xvf $ror'.tar'
gunzip -dq *.gz
gunzip -dq *.Z
#
# Xsort now works on the original fits files so don't need xrtconv
#
# produce a file with good time slots
# Work out what format data is in - us or new type rdf
set num=`echo $ror | cut -c 3-8`
set type=`cat *$num*.public_contents | grep "DATA_FORMAT" | awk '{print $3}' |cut -c 2-5`
if ( $type == "PROS" ) then
# us data format
xrthk $rawdir rootname=$ror hkpar1=iac_evr pmin1=0 pmax1=170 \
    hkpar2=iql_asp pmin2=0 pmax2=2 hkpar3=iax_evr pmin3=0 pmax3=30 \
    fname=$root'_rate_ok_1.lis' \\
else
# new type rdf
xrthk $rawdir rootname=$ror hkpar1=mv_aco pmin1=0 pmax1=170 \
    hkpar2=asp_qual pmin2=0 pmax2=2 hkpar3=xacc pmin3=0 pmax3=30 \
    fname=$root'_rate_ok_1.lis' \\
endif
#  Xsort a cube using values from defaults.txt
xsort rawdir=$rawdir rootname=$ror axes="'1 2 7'" \
    shape='c' rad=0.85 ranges="'7 5'" timrange=$root'_rate_ok_1.lis' \
    minen=$minen maxen=$maxen enbin=$enbin nxbin=$nxbin nybin=$nybin noback \
    out=$root'_cube_broad'  \\
#
#  create initial background by copying cube to bg cube
cp $root'_cube_broad.sdf' $root'_bcube_broad.sdf'
#
set imra = `hget $root'_cube_broad.more.asterix.header.axis_ra item=value'`
set imdec = `hget $root'_cube_broad.more.asterix.header.axis_dec item=value'`
set xpos = `calc exp = "'cos($dec*pi/180)*(($ra)-($imra))'"`
set ypos = `calc exp = "'($dec)-($imdec)'"`
#  save positions (in image coords) in stored.txt
if ( -e stored.txt ) rm stored.txt
echo 'xpos= '$xpos > stored.txt
echo 'ypos= '$ypos >> stored.txt
#  create point for goto to return to
redo:
#
#  mask in background cube to get bg annulus
ignore noaxsel nodatsel inp=$root'_bcube_broad' \\
crestore inp=$root'_bcube_broad' cx=0 cy=0 cr=$bgradout \\
cignore inp=$root'_bcube_broad' cx=0 cy=0 cr=$bgradin \\
#  extra ignore
source extra $root'_bcube_broad' file
#
#  create background subtracted image and bg model
xrtsub source=$root'_cube_broad' bckgnd=$root'_bcube_broad' \
    out=$root'_bscube_broad' bgmodel=$root'_bmcube_broad' \
    rootname=$ror \\
#
#  Project these down to images
project $root'_cube_broad' $root'_im_broad' axis=3 \\
project $root'_bmcube_broad' $root'_bmodel_broad' axis=3 \\
project $root'_bscube_broad' $root'_bgs_broad' axis=3 \\
#
#  search for point sources
pss inp=$root'_im_broad' bgnd=$root'_bmodel_broad' out=$root'_src_broad' \
    expert nopsfcon psf='"POLAR(XRT_PSPC,0.02)"' mask=VARPROFILE \
    aux=0.5 sigmin=4.5 map=$root'_sigmap_broad' slice='"*:*,*:*"' \
    dev=$root'_src_broad.pss' \\
#
# create ard file (contains list of regions to be ignored - e.g pt sources,
# or xrt support structure)
# xrtard is best way to do this. 
echo 'creating ardfile using xrtard'
xrtard rootname=$ror ardfile=ard_broad.dat timrange=$root'_rate_ok_1.lis' \
    srclist=$root'_src_broad' energy=0.5 pfrac=0.95 new=no mode=c \
    outside=n doplot=n \\
#
#  compare the two files with one another. If different redo pss with
#  updated background
set comp1=`wc -l ard_broad.dat | awk '{print $1}'`
set comp2=`wc -l old_ard.dat | awk '{print $1}'`
if ( $comp1 == $comp2) then
    echo 'itterations complete'
    echo 'number of itterations= '$itt
    echo ' '
else
    #  reset all data first
    quality set over inp=$root'_bcube_broad' fsel=n noaxsel nodatsel qval="good" \
    #  mask out current ard file data on background cube
    quality set over inp=$root'_bcube_broad' fsel=y noaxsel nodatsel\
       ardfile="ard_broad.dat" qval="bad"
    #  copy current ard data to old data
    cp ard_broad.dat old_ard.dat
    #  count the total no. of itterations
    @ itt=$itt + 1
    rm ard_broad.dat
    #  erase unesesary backup files
    rm *~
    #  redo the background subtraction and source serching
    goto redo
endif
#
#  now check background scaling factor
#  need to create some polar profiles first
#  only use lower energy bin(s) for this job
#  therefore need to get only lowest energy bins - need to define sensible range
set noebins=`hget inp=$root'_cube_broad.axis(3).data_array' item=dims \\`
# work with approx bottom 10%
set upper=`calc exp="'nint((($maxen-$minen)/10)+$minen)'"`
# check that this actually includes at least one bin - if not, make it so
set check=`hget inp=$root'_cube_broad.axis(3).data_array(1)' item=value \\`
set check2=`calc exp="'nint($check)'"`
if ( $check2 > $upper ) set upper=$check
#
binsubset inp=$root'_cube_broad' out=datlow axes=3 axis3=$minen':'$upper \\
binsubset inp=$root'_bmcube_broad' out=bmlow axes=3 axis3=$minen':'$upper \\
# unfortunately this may leave a cube with several energy bins - we want an image
# if not an image - make it so
set tdims=`hdir datlow.axis | grep "dimensions" | awk '{print $7}' |cut -c 2`
if ( $tdims == 3 ) then
project datlow datlow2 axis=3 \\
project bmlow bmlow2 axis=3 \\
rm datlow.sdf bmlow.sdf
mv datlow2.sdf datlow.sdf
mv bmlow2.sdf bmlow.sdf
endif
#  mask image with ard file and make polar profile
quality set over fsel=y inp=datlow axes=\"1 2\" \
    ardfile="ard_broad.dat" qval="bad" \\
#  extra ignore
source extra datlow file
#
ipolar inp=datlow out=datlow_polar xcent=0 ycent=0 \
    rbin=0.01 abin=360 \\
ipolar inp=bmlow out=bmlow_polar xcent=0 ycent=0 \
    rbin=0.01 abin=360 \\
#
#  divide polar profile by bg model profile
divide datlow_polar bmlow_polar $root'_bg_check'
ignore inp=$root'_bg_check' noaxsel nodatsel \\
restore nodatsel inp=$root'_bg_check' axis1='"'$bgradin':'$bgradout'"'
#
#  get bg scaling factor
#  first get initial mean and standard deviation.
set inmean = `statistix $root'_bg_check' | grep "Weighted mean" | awk '{print $5}'`
set indev = `statistix $root'_bg_check' | grep "Standard deviation" | awk '{print $5}'`
#  ignore points more than two standard deviations from the mean
set min = `calc exp="'$inmean - $indev'"`
set max = `calc exp="'$inmean + $indev'"`
ignore inp=$root'_bg_check' noaxsel nodatsel \\
restore inp=$root'_bg_check' noaxsel data='"'$min':'$max'"' \\
ignore inp=$root'_bg_check' nodatsel axis1='"0:'$bgradin'"' \\
ignore inp=$root'_bg_check' nodatsel axis1='"'$bgradout':0.9"' \\
set bgscale = `statistix $root'_bg_check' | grep "Weighted mean" | awk '{print $5}'`
#
set check2=`calc exp="'nint(1000*$bgscale)'"`
if ( ( $check2 > 1100 ) || ( $check2 < 900 ) ) echo 'WARNING - SCALING FACTOR MAY BE TOO BIG/SMALL'
# if value is within 2% don't do anything else mult by bgscale
if ( ( $check2 > 1020 ) || ( $check2 < 980 ) ) then
echo 'background correction = '$bgscale
multiply over $root'_bcube_broad' $bgscale
else
echo 'no background correction needed'
endif
#
#
#  final background subtraction and exposure correction of source cube
xrtsub source=$root'_cube_broad' bckgnd=$root'_bcube_broad' \
    out=$root'_bscube_broad' bgmodel=$root'_bmcube_broad' \
    rootname=$ror sdefvar=0 bdefvar=0 \\
xrtcorr inp=$root'_bscube_broad' out=$root'_ccube_broad' rtname=$ror \\
#
#  Construct spectral response by projecting the corrected cube to a
#  spectrum and then copying the response to the cube
#
project $root'_cube_broad' temp1 1
project temp1 temp2 1
xrtcorr temp2 temp3 pcorr=no offax1=0 \\
xrtresp temp3 erase=y \\
hcopy temp3.more.asterix.energy_resp $root'_ccube_broad.more.asterix.energy_resp'
rm -f temp1.sdf temp2.sdf temp3.sdf
#
#  project cubes
project $root'_cube_broad' $root'_image' 3
project $root'_ccube_broad' $root'_cimage' 3
project $root'_bmcube_broad' $root'_bimage' 3
#
cp $root'_image.sdf' temp.sdf
quality set over fsel=y inp=temp noaxsel nodatsel ardfile="ard_broad.dat" qval="bad" \\
source extra temp file
#
#  get final image of backgrnd and orig image compared
ipolar temp $root'_image_1_p' xcent=$xpos ycent=$ypos rbin=0.02 abin=360 \\
ipolar $root'_bimage' $root'_bimage_1_p' xcent=$xpos ycent=$ypos rbin=0.02 abin=360 \\
#
divide inp1=$root'_image_1_p' inp2=$root'_bimage_1_p' out=profile
#
# clean up unnecessary files
rm temp.sdf bmlow.sdf bmlow_polar.sdf datlow.sdf datlow_polar.sdf
rm old_ard.dat $root'_im_broad.sdf' $root'_bmodel_broad.sdf'
rm $root'_bgs_broad.sdf' $root'_bg_check.sdf'
#
exit


    

    
Asterix


Last modified: Fri Jun  2 14:19:31 BST 2000