#! /bin/tcsh -f

# t_spec 4.0

# Performs additional processing for spectral files (rmf etc.)
# Uses code from do_spec (v3.2) script by K.Nandra/D.Turcan

# Adaptation by Paul O'Neill

# Calls the script: t_group

# Date in final form: 17 December 2004
 

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

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


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

set obsid = $1           # observation sequence number
set datamode = $2        # dominant datamode of this sequence
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         # location of Tartarus 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


#### check for raw .pha files and files with centroid in binned detector coordinates ####

foreach inst ( sis0 sis1 gis2 gis3 )

    set file = ${obsid}_${inst}_src.pha
    if ( ! -e $WORK/$file ) then
	echo "gnt_fatal (t_spec) File $file is missing."
	exit 1
    endif

    set file = detcent_${inst}.temp     # centroids in binned detector coordinates - used by ascaarf
    if ( ! -e $WORK/$file ) then
	echo "gnt_fatal (t_spec) File $file is missing."
	exit 1
    endif

end


#### check for otime file ####

if ( ! -e $WORK/otime.temp ) then
    echo "gnt_fatal (t_spec) File otime.temp is missing."
    exit 1
else
    set otime = `head -1 otime.temp`
endif





########################################################################################
# Get the appropriate response matrices from the local caldb and rename
########################################################################################

# NOTE: caldb must be initialised so that the environment variable $CALDB is defined.  The
# script run_tartarus checks for this and exits if caldb is has not been initialised.

echo "gnt_infrm (t_spec) getting rmf files from CALDB"


#### get SIS rmf files ####

if ( $datamode == "BRIGHT2" ) then
    cp $CALDB/data/asca/sis/cpf/94nov9/s0c1g0234p40e0_1024v0_8i.rmf $WORK/${obsid}_sis0.rmf
    cp $CALDB/data/asca/sis/cpf/94nov9/s1c3g0234p40e0_1024v0_8i.rmf $WORK/${obsid}_sis1.rmf
else
    cp $CALDB/data/asca/sis/cpf/94nov9/s0c1g0234p40e1_512v0_8i.rmf $WORK/${obsid}_sis0.rmf
    cp $CALDB/data/asca/sis/cpf/94nov9/s1c3g0234p40e1_512v0_8i.rmf $WORK/${obsid}_sis1.rmf
endif


#### get GIS rmf files ####

if ( $bitfix == "y" ) then
    cp $CALDB/data/asca/gis/cpf/95mar06/gis2v4_0.rmf $WORK/${obsid}_gis2.rmf
    cp $CALDB/data/asca/gis/cpf/95mar06/gis3v4_0_128ch.rmf $WORK/${obsid}_gis3.rmf
else
    cp $CALDB/data/asca/gis/cpf/95mar06/gis2v4_0.rmf $WORK/${obsid}_gis2.rmf
    cp $CALDB/data/asca/gis/cpf/95mar06/gis3v4_0.rmf $WORK/${obsid}_gis3.rmf
endif


########################################################################################
# SIS0 and SIS1  -  alter keywords, set BAD spectral channels, and group
########################################################################################

cd $WORK

echo "gnt_infrm (t_spec) altering keywords and grouping SIS data"

foreach inst ( sis0 sis1 )

    set file = ${obsid}_${inst}_src.pha
    if ( -e ${file}.tmp ) rm -f ${file}.tmp


#### alter the keywords ####

    echo "gnt_infrm (t_spec) altering keywords for file $file"

    set grpcom = ( chkey resp ${obsid}_${inst}.rmf & chkey ancrfile ${obsid}_${inst}.arf ) 
    set grpcom = ( $grpcom & chkey backfile ${obsid}_${inst}_bgd.pha & exit ) # command for grppha

    echo "gnt_infrm (t_spec)"
    echo "gnt_infrm (t_spec) grouping file ${file}.tmp"
    echo "gnt_infrm (t_spec)"
    echo "gnt_infrm (t_spec) grppha command : $grpcom"
    echo "gnt_infrm (t_spec)"

    grppha infile=$file outfile=${file}.tmp comm="${grpcom}"     # set keywords


#### set grouping command for SIS data ####

    if ( $datamode == "BRIGHT2" ) then

	echo "gnt_infrm (t_spec)" 
	echo "gnt_infrm (t_spec) using standard BRIGHT2 grouping for $inst"
	echo "gnt_infrm (t_spec)"
	set grpcom = ( reset quality & reset grouping & bad 0-33 & bad 682-1023 )
	set grpcom = ( $grpcom & group 34 73 2 & group 74 281 4 & group 282 481 8 & group 482 553 24 )
	set grpcom = ( $grpcom & group 554 681 64 & exit )

    else     # datamode is BRIGHT

	echo "gnt_infrm (t_spec)"
	echo "gnt_infrm (t_spec) using standard BRIGHT grouping for $inst"
	echo "gnt_infrm (t_spec)"
	set grpcom = ( reset quality & reset grouping & bad 0-16 & bad 342-511 )
	set grpcom = ( $grpcom & group 42 145 2 & group 146 241 4 & group 242 277 12 & group 278 341 32 & exit ) 

    endif


    echo "gnt_infrm (t_spec)"
    echo "gnt_infrm (t_spec) grouping file ${file}.tmp"
    echo "gnt_infrm (t_spec)"
    echo "gnt_infrm (t_spec) grppha command : $grpcom"
    echo "gnt_infrm (t_spec)"

    grppha infile=${file}.tmp outfile=${obsid}_${inst}_src_std.pha clobber=yes comm="${grpcom}" # standard grouping

    $SCRIPTS/t_group ${obsid}_${inst}_src_std.pha ${obsid}_${inst}_src_20.pha 20      # group min 20
    
    if ( $status != "0" ) then
	echo "gnt_fatal (t_spec) t_group failed (sequence is $obsid)"
	echo "gnt_fatal (t_spec) exiting from t_spec"
	exit 1
    endif

    mv ${obsid}_${inst}_src.pha.tmp ${obsid}_${inst}_src.pha    # this is the ungrouped pha file

end


########################################################################################
# GIS2 and GIS3  -  alter keywords, set BAD spectral channels, and group 
########################################################################################

cd $WORK

echo "gnt_infrm (t_spec) altering keywords and grouping GIS data"

foreach inst ( gis2 gis3 )

    set file = ${obsid}_${inst}_src.pha
    if ( -e ${file}.tmp ) rm -f ${file}.tmp


#### alter the keywords ####

    set grpcom = ( chkey resp ${obsid}_${inst}.rmf & chkey ancrfile ${obsid}_${inst}.arf ) 
    set grpcom = ( $grpcom & chkey backfile ${obsid}_${inst}_bgd.pha & exit ) # command for grppha
    
    echo "gnt_infrm (t_spec)"
    echo "gnt_infrm (t_spec) altering keywords for file $file"
    echo "gnt_infrm (t_spec)"
    echo "gnt_infrm (t_spec) grppha command : $grpcom"
    echo "gnt_infrm (t_spec)"

    grppha infile=$file outfile=${file}.tmp comm="${grpcom}"     # set keywords


#### set grouping command for bitfix GIS3 data ####

    if ( ($bitfix == "y") && ($inst == "gis3") ) then    # only 128 channels for GIS3

	echo "gnt_infrm (t_spec)" 
	echo "gnt_infrm (t_spec) using bitfix (128 channels) grouping for $inst"
	echo "gnt_infrm (t_spec)"

	if ( $otime < 108000000 ) then    # for observations before mission time 108000000 (before July 1996)

	    echo "gnt_infrm (t_spec) otime is less than 108000000"
	    set grpcom = ( reset quality & reset grouping & bad 0-7 & bad 107-127 & exit )
	    	    
	else   # later observations

	    echo "gnt_infrm (t_spec) otime is greater than 108000000"
	    set grpcom = ( reset quality & reset grouping & bad 0-11 & bad 107-127 & exit )

	endif


#### set grouping command for GIS2 or non-bitfix data ####

    else     # 1024 channels if GIS2 or if bitfix is no

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

	if ( $otime < 108000000 ) then    # for observations before mission time 108000000 (before July 1996) 

	    echo "gnt_infrm (t_spec) otime is less than 108000000"
	    set grpcom = ( reset quality & reset grouping & bad 0-60 & bad 856-1023 )
	    set grpcom = ( $grpcom & group 61 102 2 & group 103 374 4 & group 375 630 8 )
	    set grpcom = ( $grpcom & group 631 726 24 & group 727 854 64 & exit )

	else   # later observations

	    echo "gnt_infrm (t_spec) otime is greater than 108000000"
 	    set grpcom = ( reset quality & reset grouping & bad 0-88 & bad 856-1023 )
	    set grpcom = ( $grpcom & group 89 102 2 & group 103 374 4 & group 375 630 8 )
	    set grpcom = ( $grpcom & group 631 726 24 & group 727 854 64 & exit )

	endif


    endif


#### do the grouping ####

    echo "gnt_infrm (t_spec)"
    echo "gnt_infrm (t_spec) grouping file ${file}.tmp"
    echo "gnt_infrm (t_spec)"
    echo "gnt_infrm (t_spec) grppha command : $grpcom"
    echo "gnt_infrm (t_spec)"

    grppha infile=${file}.tmp outfile=${obsid}_${inst}_src_std.pha clobber=yes comm="${grpcom}" # standard grouping

    $SCRIPTS/t_group ${obsid}_${inst}_src_std.pha ${obsid}_${inst}_src_20.pha 20      # group min 20
    
    if ( $status != "0" ) then
	echo "gnt_fatal (t_spec) t_group failed (sequence is $obsid)"
	echo "gnt_fatal (t_spec) exiting from t_spec"
	exit 1
    endif

    mv ${obsid}_${inst}_src.pha.tmp ${obsid}_${inst}_src.pha    # this is the ungrouped pha file

end


########################################################################################
# SIS0 - create arf files
########################################################################################

# ascaarf is used to create ancillary response files
# the object position is obtained from the detector centroid file

set srcpos = `cat $WORK/detcent_sis0.temp`

if ( "$#srcpos" != "2" ) then    # region not found
    echo "gnt_fatal (t_spec) could not get the source position for SIS0 - exiting"
    exit 1
else
    set srcpos = ( `echo "scale=3;$srcpos[1] * 4.0" | bc` `echo "scale=3;$srcpos[2] * 4.0" | bc` )
    echo "gnt_infrm (t_spec) unbinned source position (DETX DETY):  $srcpos[1] $srcpos[2]"
    set arfcmd = "point=yes simple=no Xposition=$srcpos[1] Yposition=$srcpos[2] clobber=yes"
    echo "gnt_infrm (t_spec) Creating ancillary response file for sis0"
    ascaarf ${obsid}_sis0_src.pha ${obsid}_sis0.rmf ${obsid}_sis0.arf $arfcmd
endif


########################################################################################
# SIS1 - create arf file
########################################################################################

# ascaarf is used to create ancillary response files
# the object position is obtained from the detector centroid file

set srcpos = `cat $WORK/detcent_sis1.temp`

if ( "$#srcpos" != "2" ) then    # region not found
    echo "gnt_fatal (t_spec) could not get the source position for SIS1 - exiting"
    exit 1
else
    set srcpos = ( `echo "scale=3;$srcpos[1] * 4.0" | bc` `echo "scale=3;$srcpos[2] * 4.0" | bc` )
    echo "gnt_infrm (t_spec) unbinned source position (DETX DETY):  $srcpos[1] $srcpos[2]"
    set arfcmd = "point=yes simple=no Xposition=$srcpos[1] Yposition=$srcpos[2] clobber=yes"
    echo "gnt_infrm (t_spec) Creating ancillary response file for sis1"
    ascaarf ${obsid}_sis1_src.pha ${obsid}_sis1.rmf ${obsid}_sis1.arf $arfcmd
endif


########################################################################################
# GIS2 - alter keywords and groupings in pha files, and create arf files - GIS2
########################################################################################

# ascaarf is used to create ancillary response files
# the object position is obtained from the detector centroid file

set srcpos = `cat $WORK/detcent_gis2.temp`

if ( "$#srcpos" != "2" ) then    # region not found
    echo "gnt_fatal (t_spec) could not get the source position for GIS2 - exiting"
    exit 1
else
    set srcpos = ( `echo "scale=3;$srcpos[1]" | bc` `echo "scale=3;$srcpos[2]" | bc` )
    echo "gnt_infrm (t_spec) unbinned source position (DETX DETY):  $srcpos[1] $srcpos[2]"
    set arfcmd = "point=yes simple=no Xposition=$srcpos[1] Yposition=$srcpos[2] clobber=yes"
    echo "gnt_infrm (t_spec) Creating ancillary response file for gis2"
    ascaarf ${obsid}_gis2_src.pha ${obsid}_gis2.rmf ${obsid}_gis2.arf $arfcmd
endif


########################################################################################
# GIS3 - alter keywords and groupings in pha files, and create arf files - GIS3
########################################################################################

# ascaarf is used to create ancillary response files
# the object position is obtained from the detector centroid file

set srcpos = `cat $WORK/detcent_gis3.temp`

if ( "$#srcpos" != "2" ) then    # region not found
    echo "gnt_fatal (t_spec) could not get the source position for GIS3 - exiting"
    exit 1
else
    set srcpos = ( `echo "scale=3;$srcpos[1]" | bc` `echo "scale=3;$srcpos[2]" | bc` )
    echo "gnt_infrm (t_spec) unbinned source position (DETX DETY):  $srcpos[1] $srcpos[2]"
    set arfcmd = "point=yes simple=no Xposition=$srcpos[1] Yposition=$srcpos[2] clobber=yes"
    echo "gnt_infrm (t_spec) Creating ancillary response file for gis3"
    ascaarf ${obsid}_gis3_src.pha ${obsid}_gis3.rmf ${obsid}_gis3.arf $arfcmd
endif


########################################################################################
# Create SIS rmf files
########################################################################################

echo "gnt_infrm (t_spec) creating rmf files"

foreach inst ( sis0 sis1 )

    rm ${obsid}_${inst}.rmf

    if ( -e ${obsid}_${inst}.arf ) then
	pset sisrmg phtopi=$CALDB/data/asca/sis/bcf/sisph2pi_290301.fits    # this is the latest file
	pset sisrmg ecdata=CALDB echosh=CALDB rddhis=CALDB datadir="" clobber=yes
	echo "gnt_infrm (t_spec) Creating response matrix for $inst"
	sisrmg infile=${obsid}_${inst}_src.pha arfile=${obsid}_${inst}.arf rmfile=${obsid}_${inst}.rmf
    else
	echo "gnt_fatal (t_spec) arf for $inst not found.  NOT creating new rmf file. Exiting."
	exit 1
    endif

end


########################################################################################
# Check that all files were created
########################################################################################

foreach inst ( sis0 sis1 gis2 gis3 )
    foreach ext ( .rmf .arf _src_std.pha _src_20.pha )
	set file = ${obsid}_${inst}$ext
	if ( ! -e $WORK/$file ) then
	    echo "gnt_fatal (t_spec) File $file not created."
	    exit 1
	endif
    end
end


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

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