#! /bin/tcsh -f

# t_expo 4.0

# Creates exposure maps.
# Uses code from do_expo (v3.2) script by D.Turcan

# Adaptation by Paul O'Neill

# Date in final form: 17 December 2004
 

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

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


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

set obsid = $1           # observation sequence number
set datamode = $2        # dominant datamode of this suquence
set bitfix = $3          # is there a GIS3 bitfix problem (y = bitfix detected, n = no fix detected)
set SEQ = $4             # root directory for current sequence
set SCRIPTS = $5         # directory where the standard spectrum is (in scripts dir)


########################################################################################
# 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

#### check for pha, rmf, and arf files ####


foreach inst ( sis0 sis1 gis2 gis3 )

    set file = ${obsid}_${inst}.rmf

    if ( ( $inst == sis0 || $inst == sis1 ) && ! -e $WORK/$file ) then
	echo "gnt_fatal (t_expo) File $file does not exist."
	exit 1
    endif

    foreach ext ( _src.pha .arf _src_std.pha )
	set file = ${obsid}_${inst}${ext}
	if ( ! -e $WORK/$file ) then
	    echo "gnt_fatal (t_expo) File $file does not exist."
	    exit 1
	endif
    end


end


#### check for attitude files ####

set attfiles = `ls -1 $AUX/fa*`

if ( "$#attfiles" == "0" ) then
    echo "gnt_fatal (t_expo) attitude files not found"
    exit 1
endif


#### get observation time ####

if ( ! -e $WORK/otime.temp ) then
    echo "gnt_fatal (t_expo) otime.temp not found"
else
    set otime = `cat $WORK/otime.temp`
    echo "gnt_infrm (t_expo) otime is $otime"
endif


########################################################################################
# Determine the BAD channels for the std spectra
########################################################################################

# find the channel numbers of the first and last good channels in each of the standard regrouped
# spectra. These numbers are used to create a 4-channel spectrum and also for input into
# ascaeffmap

foreach inst ( sis0 sis1 gis2 gis3 )


#### get first good channel ###

# Loops until finds a good channel

    echo "gnt_infrm (t_expo)"
    echo "gnt_infrm (t_expo) finding first good channel for $inst"

    @ row = 1

    ftabpar fitsfile="${obsid}_${inst}_src_std.pha+1" column="QUALITY" row="$row"
    set quality = `pget ftabpar value`

    while ( $quality != 0 )
	@ row ++
	ftabpar fitsfile="${obsid}_${inst}_src_std.pha+1" column="QUALITY" row="$row"
	set quality = `pget ftabpar value`
    end

    if ( $inst == sis0 ) then
	@ sis0g1 = $row - 1   # row 1 corresponds to channel 0
    endif

    if ( $inst == sis1 ) then
	@ sis1g1 = $row - 1
    endif

    if ( $inst == gis2 ) then
	@ gis2g1 = $row - 1
    endif

    if ( $inst == gis3 ) then
	@ gis3g1 = $row - 1
    endif


#### get last good channel ####

    echo "gnt_infrm (t_expo) finding last good channel for $inst"

    fkeypar fitsfile="${obsid}_${inst}_src_std.pha+1" keyword="NAXIS2"
    set row = `pget fkeypar value`

    ftabpar fitsfile="${obsid}_${inst}_src_std.pha+1" column="QUALITY" row="$row"
    set quality = `pget ftabpar value`

    while ( $quality != 0 )
	@ row --
	ftabpar fitsfile="${obsid}_${inst}_src_std.pha+1" column="QUALITY" row="$row"
	set quality = `pget ftabpar value`
    end

    if ( $inst == sis0 ) then
	@ sis0g2 = $row - 1   # row 1 corresponds to channel 0
    endif

    if ( $inst == sis1 ) then
	@ sis1g2 = $row - 1
    endif

    if ( $inst == gis2 ) then
	@ gis2g2 = $row - 1
    endif

    if ( $inst == gis3 ) then
	@ gis3g2 = $row - 1
    endif


#### show the channels ####

    if ( $inst == sis0 ) then
	echo "gnt_infrm (t_expo) SIS0 first and last good channels in spectrum: $sis0g1   $sis0g2"
    endif

    if ( $inst == sis1 ) then
	echo "gnt_infrm (t_expo) SIS1 first and last good channels in spectrum: $sis1g1   $sis1g2"
    endif

    if ( $inst == gis2 ) then
	echo "gnt_infrm (t_expo) GIS2 first and last good channels in spectrum: $gis2g1   $gis2g2"
    endif

    if ( $inst == gis3 ) then
	echo "gnt_infrm (t_expo) GIS3 first and last good channels in spectrum: $gis3g1   $gis3g2"
    endif

end


########################################################################################
# Check for chipid for SIS0 and SIS1
########################################################################################

set chipid0 = ""

foreach inst ( 0 1 )

    set chip = `fkeyprint ${obsid}_sis${inst}.evt+0 S${inst}CCDPOW | grep S${inst}CCDPOW=`
    set chip = `echo $chip[2] | awk '{print substr($1,2,4)}'`
    set chipid = ""

    foreach num ( 0 1 2 3 )
	if ( `echo $chip $num | awk '{print substr($1,$2+1,1)}'` == 1 ) \
	    set chipid = ${chipid}$num
    end

    if ( $chipid0 == "" ) then
	set chipid0 = $chipid
    else
	set chipid1 = $chipid
    endif

end


########################################################################################
# Simulate standard spectra for all instruments
########################################################################################

# Create then use an XSPEC script to simulate spectra with gamma=1.4, 1.7, 2.0. Efficiency
# maps are then created for these standard spectrum in addition to those created using
# the source spectrum.

echo "gnt_infrm (t_expo) faking standard spectra"

set gammas = ( 1.4 1.7 2.0 )

#### start loop through gamma ####

foreach gamma ( $gammas )

    echo "gnt_infrm (t_expo)"
    echo "gnt_infrm (t_expo) creating XSPEC script for gamma = $gamma"


#### remove any previous files ####

    foreach inst ( sis0 sis1 gis2 gis3 )
	if ( -e sim_std_spec_${gamma}.xcm ) rm sim_std_spec_${gamma}.xcm                   # script 
	if ( -e ${obsid}_${inst}_src.fak ) rm ${obsid}_${inst}_src.fak               # script creates these files
	if ( -e ${obsid}_${inst}_fak_${gamma}.pha ) rm ${obsid}_${inst}_fak_${gamma}.pha   # spectra are renamed to this
    end


#### create XSPEC script ####

    echo "set xs_return_result 1" >! sim_std_spec_${gamma}.xcm

    echo "data 1:1 ${obsid}_sis0_src" >> sim_std_spec_${gamma}.xcm 
    echo "data 2:2 ${obsid}_sis1_src" >> sim_std_spec_${gamma}.xcm
    echo "data 3:3 ${obsid}_gis2_src" >> sim_std_spec_${gamma}.xcm
    echo "data 4:4 ${obsid}_gis3_src" >> sim_std_spec_${gamma}.xcm

    echo "ignore bad" >> sim_std_spec_${gamma}.xcm

    echo 'model po' >> sim_std_spec_${gamma}.xcm

    echo "${gamma},0.2,-1,-1,3,3" >> sim_std_spec_${gamma}.xcm
    echo "1e2,.0002,0,0,1e4,1e4" >> sim_std_spec_${gamma}.xcm
    echo "${gamma},0.2,-1,-1,3,3" >> sim_std_spec_${gamma}.xcm
    echo "1e2,.0002,0,0,1e4,1e4" >> sim_std_spec_${gamma}.xcm
    echo "${gamma},0.2,-1,-1,3,3" >> sim_std_spec_${gamma}.xcm
    echo "1e2,.0002,0,0,1e4,1e4" >> sim_std_spec_${gamma}.xcm
    echo "${gamma},0.2,-1,-1,3,3" >> sim_std_spec_${gamma}.xcm
    echo "1e2,.0002,0,0,1e4,1e4" >> sim_std_spec_${gamma}.xcm

    echo "freeze 1-8" >> sim_std_spec_${gamma}.xcm
    
    echo "backgrnd none,none,none,none" >> sim_std_spec_${gamma}.xcm   # don't want background for faked standard spectrum

    echo "fakeit" >> sim_std_spec_${gamma}.xcm
    echo "y" >> sim_std_spec_${gamma}.xcm   # need to use Poisson stats so get a spectrum with counts (rather than counts/s)
    echo " " >> sim_std_spec_${gamma}.xcm   # all these blank lines are needed to call the fakeit command
    echo " " >> sim_std_spec_${gamma}.xcm
    echo " " >> sim_std_spec_${gamma}.xcm
    echo " " >> sim_std_spec_${gamma}.xcm
    echo " " >> sim_std_spec_${gamma}.xcm
    echo " " >> sim_std_spec_${gamma}.xcm
    echo " " >> sim_std_spec_${gamma}.xcm
    echo " " >> sim_std_spec_${gamma}.xcm
    echo " " >> sim_std_spec_${gamma}.xcm
    echo " " >> sim_std_spec_${gamma}.xcm
    echo "exit" >> sim_std_spec_${gamma}.xcm
    echo " " >> sim_std_spec_${gamma}.xcm

    if ( ! -e sim_std_spec_${gamma}.xcm ) then
	echo "gnt_fatal (t_expo) script file sim_std_spec_${gamma}.xcm was not created"
	exit 1
    endif


#### run XSPEC script ####

    echo "gnt_infrm (t_expo) running XSPEC script sim_std_spec_${gamma}.xcm"    

    xspec - sim_std_spec_${gamma}.xcm


#### rename files ####

    foreach inst ( sis0 sis1 gis2 gis3 )

	if ( ! -e ${obsid}_${inst}_src.fak ) then
	    echo "gnt_fatal (t_expo) file ${obsid}_${inst}_src.fak was not created"
	    exit 1
	endif

	echo "gnt_infrm (t_expo)"

	echo "gnt_infrm (t_expo) copying ${obsid}_${inst}_src.pha to ${obsid}_${inst}_fak_${gamma}.pha"
	cp ${obsid}_${inst}_src.pha ${obsid}_${inst}_fak_${gamma}.pha         # need to do this because ascaeffmap needs a pha file with a WMAP

	echo "gnt_infrm (t_expo) removing COUNTS from ${obsid}_${inst}_fak_${gamma}.pha"
	fdelcol infile="${obsid}_${inst}_fak_${gamma}.pha+1" colname="COUNTS" confirm="no" proceed="yes"

	echo "gnt_infrm (t_expo) copying in fake COUNTS from ${obsid}_${inst}_src.fak"
	faddcol infile="${obsid}_${inst}_fak_${gamma}.pha+1" colfile="${obsid}_${inst}_src.fak+1" colname="COUNTS"

    end


#### end loop through gamma ####

end


echo "gnt_infrm (t_expo) created all fake spectra"


########################################################################################
# SIS instruments setup - group into 4 channel spectrum and make efficiency map
########################################################################################

# use the source spectrum then fasked spectra

# group channels to get a 4-channel spectrum and
# create efficiency map using ascaeffmap with these temporary spectral files

foreach inst ( s0 s1 )


#### remove previous files ####

    if ( -e ${obsid}_si${inst}_4_ch.pha ) rm ${obsid}_si${inst}_4_ch.pha
    if ( -e ${obsid}_si${inst}.effmap) rm ${obsid}_si${inst}.effmap

    foreach gamma ( $gammas )
	if ( -e ${obsid}_si${inst}_fak_${gamma}_4_ch.pha ) rm ${obsid}_si${inst}_fak_${gamma}_4_ch.pha
	if ( -e ${obsid}_si${inst}_${gamma}.effmap ) rm ${obsid}_si${inst}_${gamma}.effmap
    end


#### set some parameters ####

    if ( $inst == "s0" ) then 
	set c = c1
	set chipid = $chipid0

    endif

    if ( $inst == "s1" ) then
	set c = c3
	set chipid = $chipid1
    endif

    echo "gnt_infrm (t_expo) Creating 4-channel grouped spectra and efficiency maps for si$inst"


#### for BRIGHT mode data ####

    if ( $datamode == "BRIGHT" ) then

	echo "gnt_infrm (t_expo)"
	echo "gnt_infrm (t_expo) grouping source spectrum for si$inst"

	if ( $inst == s0 ) then
	    @ lastbad = $sis0g1 - 1       # last bad channel at start of spectrum
	    @ firstbad = $sis0g2 + 1      # first bad channel at end of spectrum
	    @ firstgood = $sis0g1
	    @ lastgood = $sis0g2
	    @ chans1 = (67 - ${sis0g1}) + 1   # number of channels in first re-grouped channel
	    @ chans2 = (${sis0g2} - 241) + 1  # number of channels in last re-grouped channel 
	    @ phalo = $firstgood + 1      # lowest regrouped channel number - for ascaeffmap
	    @ phahi = $phalo + 3          # highest regrouped channel number - for ascaeffmap
	else
	    @ lastbad = $sis1g1 - 1       # last bad channel at start of spectrum
	    @ firstbad = $sis1g2 + 1      # first bad channel at end of spectrum
	    @ firstgood = $sis0g1
	    @ lastgood = $sis0g2
	    @ chans1 = (67 - ${sis1g1}) + 1   # number of channels in first re-grouped channel
	    @ chans2 = (${sis1g2} - 241) + 1  # number of channels in last re-grouped channel 
	    @ phalo = $firstgood + 1      # lowest regrouped channel number - for ascaeffmap
	    @ phahi = $phalo + 3          # highest regrouped channel number - for ascaeffmap
	endif

	if ( ($chans1 <= 0) || ($chans2 <= 0) ) then
	    echo "gnt_fatal (t_expo) problem with grouping spectra"
	    exit 1
	endif

	set grpcom = ( reset quality & reset grouping & bad 0-${lastbad} & bad ${firstbad}-511 & group ${firstgood} 67 ${chans1} )
	set grpcom = ( $grpcom & group 68 171 104 & group 172 240 69 & group 241 ${lastgood} ${chans2} & exit )  

	echo "gnt_infrm (t_expo) grppha command: $grpcom" 
	grppha infile=${obsid}_si${inst}_src.pha outfile=${obsid}_si${inst}_4_ch.pha comm="${grpcom}"

	echo "gnt_infrm (t_expo) Creating efficiency map for si$inst"

	ascaeffmap phafile=${obsid}_si${inst}_4_ch.pha chip=$chipid pha_lo=$phalo pha_hi=$phahi \
	    rmffile=$WORK/${obsid}_si${inst}.rmf outfile=$WORK/${obsid}_si${inst}.effmap

	
	foreach gamma ( $gammas )      # now use the faked spectra

	    echo "gnt_infrm (t_expo)"
	    echo "gnt_infrm (t_expo) grouping spectrum for gamma = $gamma for si$inst"

	    set grpcom = ( reset quality & reset grouping & bad 0-${lastbad} & bad ${firstbad}-511 & group ${firstgood} 67 ${chans1} )
	    set grpcom = ( $grpcom & group 68 171 104 & group 172 240 69 & group 241 ${lastgood} ${chans2} & exit )  

	    echo "gnt_infrm (t_expo) grppha command: $grpcom" 
	    grppha infile=${obsid}_si${inst}_fak_${gamma}.pha outfile=${obsid}_si${inst}_fak_${gamma}_4_ch.pha comm="${grpcom}"

	    echo "gnt_infrm (t_expo) Creating efficiency map for gamma = $gamma for si$inst"

	    ascaeffmap phafile=${obsid}_si${inst}_fak_${gamma}_4_ch.pha chip=$chipid pha_lo=$phalo pha_hi=$phahi \
		rmffile=$WORK/${obsid}_si${inst}.rmf outfile=$WORK/${obsid}_si${inst}_${gamma}.effmap	    

	end


#### for BRIGHT2 mode data ####

    else

    	echo "gnt_infrm (t_expo)"
	echo "gnt_infrm (t_expo) grouping source spectrum for si$inst"

	if ( $inst == s0 ) then
	    @ lastbad = $sis0g1 - 1       # last bad channel at start of spectrum
	    @ firstbad = $sis0g2 + 1      # first bad channel at end of spectrum
	    @ firstgood = $sis0g1
	    @ lastgood = $sis0g2
	    @ chans1 = (136 - ${sis0g1}) + 1   # number of channels in first re-grouped channel
	    @ chans2 = (${sis0g2} - 478) + 1  # number of channels in last re-grouped channel 
	    @ phalo = $firstgood + 1      # lowest regrouped channel number - for ascaeffmap
	    @ phahi = $phalo + 3          # highest regrouped channel number - for ascaeffmap
	else
	    @ lastbad = $sis1g1 - 1       # last bad channel at start of spectrum
	    @ firstbad = $sis1g2 + 1      # first bad channel at end of spectrum
	    @ firstgood = $sis0g1
	    @ lastgood = $sis0g2
	    @ chans1 = (136 - ${sis1g1}) + 1   # number of channels in first re-grouped channel
	    @ chans2 = (${sis1g2} - 478) + 1  # number of channels in last re-grouped channel 
	    @ phalo = $firstgood + 1      # lowest regrouped channel number - for ascaeffmap
	    @ phahi = $phalo + 3          # highest regrouped channel number - for ascaeffmap
	endif

	if ( ($chans1 <= 0) || ($chans2 <= 0) ) then
	    echo "gnt_fatal (t_expo) problem with grouping spectra"
	    exit 1
	endif

	set grpcom = ( reset quality & reset grouping & bad 0-${lastbad} & bad ${firstbad}-1023 & group ${firstgood} 136 ${chans1} )
	set grpcom = ( $grpcom & group 137 340 204 & group 341 477 137 & group 478 ${lastgood} ${chans2} & exit )
    
	echo "gnt_infrm (t_expo) grppha command: $grpcom" 
	grppha infile=${obsid}_si${inst}_src.pha outfile=${obsid}_si${inst}_4_ch.pha comm="${grpcom}"

	echo "gnt_infrm (t_expo) Creating efficiency map for si$inst"

	ascaeffmap phafile=${obsid}_si${inst}_4_ch.pha chip=$chipid pha_lo=${phalo} pha_hi=${phahi} \
	    rmffile=$WORK/${obsid}_si${inst}.rmf outfile=$WORK/${obsid}_si${inst}.effmap 

	foreach gamma ( $gammas )      # now use the faked spectra

	    echo "gnt_infrm (t_expo)"
	    echo "gnt_infrm (t_expo) grouping spectrum for gamma = $gamma for si$inst"

	    set grpcom = ( reset quality & reset grouping & bad 0-${lastbad} & bad ${firstbad}-1023 & group ${firstgood} 136 ${chans1} )
	    set grpcom = ( $grpcom & group 137 340 204 & group 341 477 137 & group 478 ${lastgood} ${chans2} & exit )

	    echo "gnt_infrm (t_expo) grppha command: $grpcom" 
	    grppha infile=${obsid}_si${inst}_fak_${gamma}.pha outfile=${obsid}_si${inst}_fak_${gamma}_4_ch.pha comm="${grpcom}"

	    echo "gnt_infrm (t_expo) Creating efficiency map for gamma = $gamma for  si$inst"

	    ascaeffmap phafile=${obsid}_si${inst}_fak_${gamma}_4_ch.pha chip=$chipid pha_lo=${phalo} pha_hi=${phahi} \
		rmffile=$WORK/${obsid}_si${inst}.rmf outfile=$WORK/${obsid}_si${inst}_${gamma}.effmap 

	end


   endif


end


########################################################################################
# GIS instruments setup - group into 4 channel spectrum and make efficiency map
########################################################################################

# source spectrum


foreach inst ( gis2 gis3 )


#### remove previous files ####

    if ( -e ${obsid}_${inst}_4_ch.pha ) rm ${obsid}_${inst}_4_ch.pha
    if ( -e ${obsid}_${inst}.effmap) rm ${obsid}_${inst}.effmap

    foreach gamma ( $gammas )
	if ( -e ${obsid}_${inst}_fak_${gamma}_4_ch.pha ) rm ${obsid}_${inst}_fak_${gamma}_4_ch.pha
	if ( -e ${obsid}_${inst}_${gamma}.effmap) rm ${obsid}_${inst}_${gamma}.effmap
    end

    echo "gnt_infrm (t_expo)"
    echo "gnt_infrm (t_expo)"
    echo "gnt_infrm (t_expo) Creating 4-channel grouped spectrum and efficieny maps for $inst"


#### GIS3 data with bitfix problem ####

    if ( $bitfix == "y" && $inst == "gis3" ) then

	echo "gnt_infrm (t_expo) bitfix detected"

	echo "gnt_infrm (t_expo)"
	echo "gnt_infrm (t_expo) grouping source spectrum for ${inst}"

	@ lastbad = $gis3g1 - 1       # last bad channel at start of spectrum
	@ firstbad = $gis3g2 + 1      # first bad channel at end of spectrum
	@ firstgood = $gis3g1
	@ lastgood = $gis3g2
	@ chans1 = (22 - ${gis3g1}) + 1   # number of channels in first re-grouped channel
	@ chans2 = (${gis3g2} - 76) + 1  # number of channels in last re-grouped channel 
	@ phalo = $firstgood + 1      # lowest regrouped channel number - for ascaeffmap
	@ phahi = $phalo + 3          # highest regrouped channel number - for ascaeffmap

	if ( ($chans1 <= 0) || ($chans2 <= 0) ) then
	    echo "gnt_fatal (t_expo) problem with grouping spectra"
	    exit 1
	endif

	set grpcom = ( reset quality & reset grouping & bad 0-${lastbad} & bad ${firstbad}-127 & group ${firstgood} 22 ${chans1} )
	set grpcom = ( $grpcom & group 23 54 32 & group 55 75 21 & group 76 ${lastgood} ${chans2} & exit )

	echo "gnt_infrm (t_expo) grppha command: $grpcom" 
	grppha infile=${obsid}_${inst}_src.pha outfile=${obsid}_${inst}_4_ch.pha comm="${grpcom}" 

	echo "gnt_infrm (t_expo) Creating efficiency map for source spectrum for $inst"
      
	ascaeffmap phafile=${obsid}_${inst}_4_ch.pha pha_lo=${phalo} pha_hi=${phahi} \
		rmffile=$WORK/${obsid}_${inst}.rmf outfile=$WORK/${obsid}_${inst}.effmap        	    


	foreach gamma ( $gammas )

	    echo "gnt_infrm (t_expo)"
	    echo "gnt_infrm (t_expo) grouping spectum for gamma = $gamma for $inst"

	    set grpcom = ( reset quality & reset grouping & bad 0-${lastbad} & bad ${firstbad}-127 & group ${firstgood} 22 ${chans1} )
	    set grpcom = ( $grpcom & group 23 54 32 & group 55 75 21 & group 76 ${lastgood} ${chans2} & exit )

	    echo "gnt_infrm (t_expo) grppha command: $grpcom" 
	    grppha infile=${obsid}_${inst}_fak_${gamma}.pha outfile=${obsid}_${inst}_fak_${gamma}_4_ch.pha comm="${grpcom}" 

	    echo "gnt_infrm (t_expo) Creating efficiency map for gamma = $gamma  for $inst"
      
	    ascaeffmap phafile=${obsid}_${inst}_fak_${gamma}_4_ch.pha pha_lo=${phalo} pha_hi=${phahi} \
		    rmffile=$WORK/${obsid}_${inst}.rmf outfile=$WORK/${obsid}_${inst}_${gamma}.effmap       
		
	end


#### GIS2 data and non-bitfix GIS3 data ####

    else

	echo "gnt_infrm (t_spec) using standard (1024 channels) grouping for $inst"

	echo "gnt_infrm (t_expo)"
	echo "gnt_infrm (t_expo) grouping source spectrum for ${inst}"

	if ( $inst == gis2 ) then
	    @ lastbad = $gis2g1 - 1       # last bad channel at start of spectrum
	    @ firstbad = $gis2g2 + 1      # first bad channel at end of spectrum
	    @ firstgood = $gis2g1
	    @ lastgood = $gis2g2
	    @ chans1 = (170 - ${gis2g1}) + 1   # number of channels in first re-grouped channel
	    @ chans2 = (${gis2g2} - 595) + 1  # number of channels in last re-grouped channel 
	    @ phalo = $firstgood + 1      # lowest regrouped channel number - for ascaeffmap
	    @ phahi = $phalo + 3          # highest regrouped channel number - for ascaeffmap
	else
	    @ lastbad = $gis3g1 - 1       # last bad channel at start of spectrum
	    @ firstbad = $gis3g2 + 1      # first bad channel at end of spectrum
	    @ firstgood = $gis3g1
	    @ lastgood = $gis3g2
	    @ chans1 = (170 - ${gis3g1}) + 1   # number of channels in first re-grouped channel
	    @ chans2 = (${gis3g2} - 595) + 1  # number of channels in last re-grouped channel 
	    @ phalo = $firstgood + 1      # lowest regrouped channel number - for ascaeffmap
	    @ phahi = $phalo + 3          # highest regrouped channel number - for ascaeffmap
	endif

	if ( ($chans1 <= 0) || ($chans2 <= 0) ) then
	    echo "gnt_fatal (t_expo) problem with grouping spectra"
	    exit 1
	endif

	set grpcom = ( reset quality & reset grouping & bad 0-${lastbad} & bad ${firstbad}-1023 & group ${firstgood} 170 ${chans1} )
	set grpcom = ( $grpcom & group 171 425 255 & group 426 594 169 & group 595 ${lastgood} ${chans2} & exit )

	echo "gnt_infrm (t_expo) grppha command: $grpcom" 
	grppha infile=${obsid}_${inst}_src.pha outfile=${obsid}_${inst}_4_ch.pha comm="${grpcom}" 
      
	echo "gnt_infrm (t_expo) Creating efficiency map for source spectrum for $inst" 

	ascaeffmap phafile=${obsid}_${inst}_4_ch.pha pha_lo=${phalo} pha_hi=${phahi} \
	    rmffile=$WORK/${obsid}_${inst}.rmf outfile=$WORK/${obsid}_${inst}.effmap 
   

	foreach gamma ( $gammas )

	    echo "gnt_infrm (t_expo)"
	    echo "gnt_infrm (t_expo) grouping spectrum for gamma = $gamma for ${inst}"
	
	    set grpcom = ( reset quality & reset grouping & bad 0-${lastbad} & bad ${firstbad}-1023 & group ${firstgood} 170 ${chans1} )
	    set grpcom = ( $grpcom & group 171 425 255 & group 426 594 169 & group 595 ${lastgood} ${chans2} & exit )

	    echo "gnt_infrm (t_expo) grppha command: $grpcom" 
	    grppha infile=${obsid}_${inst}_fak_${gamma}.pha outfile=${obsid}_${inst}_fak_${gamma}_4_ch.pha comm="${grpcom}" 
      
	    echo "gnt_infrm (t_expo) Creating efficiency map for gamma = $gamma for $inst" 

	    ascaeffmap phafile=${obsid}_${inst}_fak_${gamma}_4_ch.pha pha_lo=${phalo} pha_hi=${phahi} \
		rmffile=$WORK/${obsid}_${inst}.rmf outfile=$WORK/${obsid}_${inst}_${gamma}.effmap 
       
	end


    endif

end


########################################################################################
# Check .effmap and .evt files are present
########################################################################################

foreach inst ( sis0 sis1 gis2 gis3 )

    if ( ! -e $WORK/${obsid}_${inst}.effmap ) then
	echo "gnt_fatal (t_expo) .effmap file for $inst does not exist"
	exit 1
    endif

    if ( ! -e $WORK/${obsid}_${inst}.evt ) then
	echo "gnt_fatal (t_expo) events file for $inst does not exist"
	exit 1
    endif

    foreach gamma ( $gammas )
	if ( ! -e $WORK/${obsid}_${inst}_${gamma}.effmap ) then
	    echo "gnt_fatal (t_expo) .effmap file for gamma = $gamma for $inst does not exist"
	    exit 1
	endif
    end

end


########################################################################################
# Set the attitude file
########################################################################################

cd $AUX

set attfiles = `ls -1 fa*`

if ( $#attfiles > 1 ) then

    echo "gnt_infrm (t_expo) Fmerging $#attfiles attitude files into 1"

    fmerge "$attfiles" $attfiles[1] "TIME,QPARAM,SIGMA,SENSOR" lastkey="MTIME1 MTIME1-1" clobber=yes

    foreach attfile ( $attfiles ) 
	if ( $attfile != $attfiles[1] ) rm -f $attfile
    end
   
    set attfile = $AUX/$attfiles[1]    # used the merged file and concat the directory path

else if ( $#attfiles == 1 ) then

    set attfile = $AUX/$attfiles[1]   # concat path

else

    echo "gnt_fatal (t_expo) Attitude file does not exist."
    exit 1

endif


########################################################################################
# Create exposure maps for all instruments
########################################################################################

cd $WORK

foreach inst ( sis0 sis1 gis2 gis3 )

    echo "gnt_infrm (t_expo)"
    echo "gnt_infrm (t_expo) Creating exposure map for source spectrum for $inst"
    echo "gnt_infrm (t_expo) instfile=$WORK/${obsid}_${inst}.effmap"
    echo "gnt_infrm (t_expo) Attitude file is $attfile"

    ascaexpo datafile=$WORK/${obsid}_${inst}.evt calfile=CALDB attitude=$attfile \
	ratefile=none instfile=$WORK/${obsid}_${inst}.effmap expofile=$WORK/${obsid}_${inst}.expo \
	imapfile=none skyfile=$WORK/${obsid}_${inst}_sky_raw.img attstep=15 rebin=8 binoff=0 imath=mul


    foreach gamma ( $gammas )

	echo "gnt_infrm (t_expo)"
	echo "gnt_infrm (t_expo) Creating exposure map for gamma = $gamma for $inst"
	echo "gnt_infrm (t_expo) instfile=$WORK/${obsid}_${inst}_${gamma}.effmap"
	echo "gnt_infrm (t_expo) Attitude file is $attfile"

	ascaexpo datafile=$WORK/${obsid}_${inst}.evt calfile=CALDB attitude=$attfile \
	    ratefile=none instfile=$WORK/${obsid}_${inst}_${gamma}.effmap expofile=$WORK/${obsid}_${inst}_${gamma}.expo \
	    imapfile=none skyfile=$WORK/${obsid}_${inst}_sky_raw_${gamma}.img attstep=15 rebin=8 binoff=0 imath=mul

    end


end


########################################################################################
# Check if all raw.expo and sky-image files were created
########################################################################################

set expostatus = "GOOD"

foreach inst ( sis0 sis1 gis2 gis3 )

    if ( ! -e $WORK/${obsid}_${inst}.expo ) then
	echo "gnt_fatal (t_expo) failed to create exposure for $inst."
	set expostatus = "ERROR"
    endif
    
    if ( ! -e $WORK/${obsid}_${inst}_sky_raw.img ) then
	echo "gnt_fatal (t_expo) Failed to create raw sky-image file for $inst."
	set expostatus = "ERROR"
    endif


    foreach gamma ( $gammas )

	if ( ! -e $WORK/${obsid}_${inst}_${gamma}.expo ) then
	    echo "gnt_fatal (t_expo) failed to create exposure for gamma = $gamma for  $inst."
	    set expostatus = "ERROR"
	endif
	
	if ( ! -e $WORK/${obsid}_${inst}_sky_raw_${gamma}.img ) then
	    echo "gnt_fatal (t_expo) Failed to create raw sky-image file for gamma - $gamma for $inst."
	    set expostatus = "ERROR"
	endif

    end

end


if ( $expostatus == "ERROR" ) then
    echo "gnt_fatal (t_expo) t_expo not successful - exiting"
    exit 1
endif


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

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