#! /bin/tcsh -f

# t_cols 4.0

# Determines mean counting rates for the sequences, in various bands.

# Fundamentally different from original do_cols script.

# Paul O'Neill

# Date in final form: 17 December 2004


########################################################################################
# Welcome 
########################################################################################

echo "gnt_infrm (t_cols)"
echo "gnt_infrm (t_cols) ++++++++++++++++++++++++++++++++++++++++++++++++++++++"
echo "gnt_infrm (t_cols)  t_cols 4.0  STARTED "`date`
echo "gnt_infrm (t_cols) ++++++++++++++++++++++++++++++++++++++++++++++++++++++"
echo "gnt_infrm (t_cols)"


########################################################################################
# Get parameters
########################################################################################

set obsid = $1           # observation sequence number
set SEQ = $2             # root directory for current sequence
set SCRIPTS = $3         # directory containing scripts


########################################################################################
# Set directory variables
########################################################################################

# NOTE: aux, screened, unscreened, and work directories MUST be sub-directories
# of the sequence directory (SEQ). This structure is hard wired.

set AUX = $SEQ/aux            # contains house keeping data
set SCR = $SEQ/screened       # contains screened data
set UNSCR = $SEQ/unscreened   # contains unscreened data
set WORK = $SEQ/work          # a work area


########################################################################################
# Check for necessary files
########################################################################################

cd $WORK


#### look for the spectral files ####

foreach inst ( sis0 sis1 gis2 gis3 )
foreach reg ( src_20 bgd )
    set file = ${obsid}_${inst}_${reg}.pha
    if ( ! -e $WORK/$file ) then
        echo "gnt_fatal (t_cols) File $file does not exist."
	exit 1
    endif
end
end


#### look for the required sis event files ####

foreach inst ( sis0 sis1 )
foreach reg ( src bgd )
foreach ext ( .evt _soft.evt _hard.evt _full2.evt _soft2.evt )
    set file = ${obsid}_${inst}_${reg}${ext}
    if ( ! -e $WORK/$file ) then
        echo "gnt_fatal (t_cols) File $file does not exist."
	exit 1
    endif
end
end
end


#### look for the required gis event files ####

foreach inst ( gis2 gis3 )
foreach reg ( src bgd )
foreach ext ( .evt _hard.evt )
    set file = ${obsid}_${inst}_${reg}${ext}
    if ( ! -e $WORK/$file ) then
        echo "gnt_fatal (t_cols) File $file does not exist."
	exit 1
    endif
end
end
end


#### remove old file ####

if ( -e counts_rate.temp ) rm counts_rate.temp
touch counts_rate.temp


########################################################################################
# Determine the scaling factor
########################################################################################

echo "gnt_infrm (t_cols) Determining scaling factor"


#### read the BACKSCAL keyword from the spectral files ####

echo "gnt_infrm (t_cols) getting BACKSCAL values"

fkeypar fitsfile=${obsid}_sis0_src_20.pha+1 keyword="BACKSCAL"
set src0_scal = `pget fkeypar value`

fkeypar fitsfile=${obsid}_sis1_src_20.pha+1 keyword="BACKSCAL"
set src1_scal = `pget fkeypar value`

fkeypar fitsfile=${obsid}_gis2_src_20.pha+1 keyword="BACKSCAL"
set src2_scal = `pget fkeypar value`

fkeypar fitsfile=${obsid}_gis3_src_20.pha+1 keyword="BACKSCAL"
set src3_scal = `pget fkeypar value`


fkeypar fitsfile=${obsid}_sis0_bgd.pha+1 keyword="BACKSCAL"
set bgd0_scal = `pget fkeypar value`

fkeypar fitsfile=${obsid}_sis1_bgd.pha+1 keyword="BACKSCAL"
set bgd1_scal = `pget fkeypar value`

fkeypar fitsfile=${obsid}_gis2_bgd.pha+1 keyword="BACKSCAL"
set bgd2_scal = `pget fkeypar value`

fkeypar fitsfile=${obsid}_gis3_bgd.pha+1 keyword="BACKSCAL"
set bgd3_scal = `pget fkeypar value`


#### calculate ratios between source and background BACKSCAL ####

echo "gnt_infrm (t_cols) Calculating ratios"

set r0 = `echo $src0_scal $bgd0_scal | awk '{printf "%6.5e\n", ($1/$2)}'`
set r1 = `echo $src1_scal $bgd1_scal | awk '{printf "%6.5e\n", ($1/$2)}'`
set r2 = `echo $src2_scal $bgd2_scal | awk '{printf "%6.5e\n", ($1/$2)}'`
set r3 = `echo $src3_scal $bgd3_scal | awk '{printf "%6.5e\n", ($1/$2)}'`

echo "gnt_infrm (t_cols) SIS0 src and bgd BACKSCAL and ratio : $src0_scal    $bgd0_scal    $r0"
echo "gnt_infrm (t_cols) SIS1 src and bgd BACKSCAL and ratio : $src1_scal    $bgd1_scal    $r1"
echo "gnt_infrm (t_cols) GIS2 src and bgd BACKSCAL and ratio : $src2_scal    $bgd2_scal    $r2"
echo "gnt_infrm (t_cols) GIS3 src and bgd BACKSCAL and ratio : $src3_scal    $bgd3_scal    $r3"


########################################################################################
# Determine SIS counting rates
########################################################################################

# The number of events in the relevent event files is determined and the scaled
# background is subtracted. The counting rate is determined using the exposure time.

foreach inst ( sis0 sis1 )
foreach ext ( .evt _soft.evt _hard.evt _full2.evt _soft2.evt )       # bands to use

    if ( $inst == sis0 ) set ratio = $r0
    if ( $inst == sis1 ) set ratio = $r1

    fkeypar fitsfile=${obsid}_${inst}_src${ext}+1 keyword=NAXIS2    # get number of src region events
    set sb_evt = `pget fkeypar value`
    echo "gnt_infrm (t_cols) $inst $ext src region counts: $sb_evt"

    fkeypar fitsfile=${obsid}_${inst}_bgd${ext}+1 keyword=NAXIS2    # get number of bgd region events
    set b_totevt = `pget fkeypar value`
    echo "gnt_infrm (t_cols) $inst $ext bgd region counts: $b_totevt"

    set b_evt = `echo $b_totevt $ratio | awk '{printf "%10.9e\n", ($1*$2)}'`          # scale by BACKSCAL ratio
    echo "gnt_infrm (t_cols) $inst $ext scaled bgd region counts: $b_evt"

    set s_evt = `echo $sb_evt $b_evt | awk '{printf "%10.9e\n", ($1-$2)}'`             # background subtracted
    echo "gnt_infrm (t_cols) $inst $ext background-subtracted counts: $s_evt"

    fkeypar fitsfile=${obsid}_${inst}_src${ext}+2 keyword=EXPOSURE        # total exposure time
    set exposure = `pget fkeypar value`
    echo "gnt_infrm (t_cols) $inst $ext exposure time (s): $exposure"

    set rate = `echo $s_evt $exposure | awk '{printf "%10.9e\n", ($1/$2)}'`              # counting rate
    echo "gnt_infrm (t_cols) $inst $ext source only counting rate (counts/s): $rate"

    set b_evt = `echo $b_evt | awk '{printf "%1.3f\n", $1}'`    # 3 decimal places
    set s_evt = `echo $s_evt | awk '{printf "%1.3f\n", $1}'`    # 3 decimal places
    set rate = `echo $rate | awk '{printf "%3.2e\n", $1}'`      # 3 significant figures

    echo "gnt_infrm (t_cols) writing: $sb_evt $b_totevt $b_evt $s_evt $rate"
    echo "$sb_evt $b_totevt $b_evt $s_evt $rate" >> counts_rate.temp

    echo "gnt_infrm (t_cols)"

end
end


########################################################################################
# Determine GIS counting rates
########################################################################################

# The number of events in the relevent event files is determined and the scaled
# background is subtracted. The counting rate is determined using the exposure time.

foreach inst ( gis2 gis3 )
foreach ext ( .evt _hard.evt )

    if ( $inst == gis2 ) set ratio = $r2
    if ( $inst == gis3 ) set ratio = $r3

    fkeypar fitsfile=${obsid}_${inst}_src${ext}+1 keyword=NAXIS2    # get number of src region events
    set sb_evt = `pget fkeypar value`
    echo "gnt_infrm (t_cols) $inst $ext src region counts: $sb_evt"

    fkeypar fitsfile=${obsid}_${inst}_bgd${ext}+1 keyword=NAXIS2    # get number of bgd region events
    set b_totevt = `pget fkeypar value`
    echo "gnt_infrm (t_cols) $inst $ext bgd region counts: $b_totevt"

    set b_evt = `echo $b_totevt $ratio | awk '{printf "%10.9e\n", ($1*$2)}'`          # scale by BACKSCAL ratio
    echo "gnt_infrm (t_cols) $inst $ext scaled bgd region counts: $b_evt"

    set s_evt = `echo $sb_evt $b_evt | awk '{printf "%10.9e\n", ($1-$2)}'`             # background subtracted
    echo "gnt_infrm (t_cols) $inst $ext background-subtracted counts: $s_evt"

    fkeypar fitsfile=${obsid}_${inst}_src${ext}+2 keyword=EXPOSURE        # total exposure time
    set exposure = `pget fkeypar value`
    echo "gnt_infrm (t_cols) $inst $ext exposure time (s): $exposure"

    set rate = `echo $s_evt $exposure | awk '{printf "%10.9e\n", ($1/$2)}'`              # counting rate
    echo "gnt_infrm (t_cols) $inst $ext source only counting rate (counts/s): $rate"

    set b_evt = `echo $b_evt | awk '{printf "%1.3f\n", $1}'`    # 3 decimal places
    set s_evt = `echo $s_evt | awk '{printf "%1.3f\n", $1}'`    # 3 decimal places
    set rate = `echo $rate | awk '{printf "%3.2e\n", $1}'`      # 3 significant figures

    echo "gnt_infrm (t_cols) writing: $sb_evt $b_totevt $b_evt $s_evt $rate"
    echo "$sb_evt $b_totevt $b_evt $s_evt $rate" >> counts_rate.temp

    echo "gnt_infrm (t_cols)"

end
end


########################################################################################
# Check if all files were created
########################################################################################

cd $WORK


#### check if file exists ####

if ( ! -e counts_rate.temp ) then
    echo "gnt_fatal (t_cols) File counts_rate.temp was not created."
    exit 1
endif


#### check if correct number of lines ####

set numline = `wc counts_rate.temp`

if ( $numline[1] != 14 ) then
    echo "gnt_fatal (t_cols) File counts_rate.temp has the wrong number of lines"
    exit 1
endif


########################################################################################
# All done 
########################################################################################

echo "gnt_infrm (t_cols)"
echo "gnt_infrm (t_cols) ++++++++++++++++++++++++++++++++++++++++++++++++++++++"
echo "gnt_infrm (t_cols)  t_cols 4.0  FINISHED "`date`
echo "gnt_infrm (t_cols) ++++++++++++++++++++++++++++++++++++++++++++++++++++++"
echo "gnt_infrm (t_cols)"