#! /bin/tcsh -f

# t_fits 4.0

# Performs power-law spectral fitting using Xspec. 

# Uses code from do_fits (v3.2) script by D.Turcan

# Calls the scripts:

# Uses the scripts: ezyfit_v4.tcl t_fits_pl_fullf.xcm t_fits_pl_softf.xcm
#                   t_fits_pl_hardf.xcm t_fits_pl_nhfullf.xcm
#                   t_fits_pl_nhsoftf.xcm t_fits_pl_nhhardf.xcm

# Adaptation by Paul O'Neill

# Date in final form: 17 December 2004


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

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


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

set obsid = $1           # observation sequence number
set SEQ = $2             # root directory for current sequence
set XCO = $3             # directory containing Xspec scripts
set SCRIPTS = $4         # directory containing 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 if object was detected  - start of "if object detected" statement
########################################################################################

cd $WORK

if ( -e sds.temp ) then
    set s = `head -1 sds.temp`
    if ( $#s == "0" ) then
	echo "gnt_fatal (t_fits) the sds.temp file is empty"
	exit 1
    endif
else
    echo "gnt_fatal (t_fits) the sds.temp file does not exist"
    exit 1
endif

set detect = `echo $s | awk '{print ($1 >= 5.0)}'`  # = 1 if $s >= 5.0

if ( $detect == 1 ) then

    echo "gnt_infrm (t_fits) object detected with significance $s"

# now do the fitting


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

#### check for source spectral files ####

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


#### check for response matrices and ancillary response files ####

foreach inst ( sis0 sis1 gis2 gis3 )
    foreach ext ( arf rmf )
	set file = ${obsid}_${inst}.$ext
	if ( ! -e $WORK/$file || `test ! -s $WORK/$file && echo no` == no ) then
	    echo "gnt_fatal (t_fits) file $file does not exist or is empty."
	    exit 1
	endif
    end
end


########################################################################################
# Run Xspec
########################################################################################


#### remove previous fits ####

if ( -e $WORK/plfit.log ) rm -r $WORK/plfit.log 

foreach band ( fullf softf hardf nhfullf nhsoftf nhhardf )

    if ( -e ${obsid}_webfit_${band}.xcm ) rm -f ${obsid}_webfit_${band}.xcm
    if ( -e fit_result_pl_${band}.dat ) rm -f fit_result_pl_${band}.dat
    if ( -e fit_result_flux_${band}.dat ) rm -f fit_result_flux_${band}.dat
    if ( -e plfit_${band}.log ) rm -f plfit.log
    if ( -e web_pha_${band}.gif ) rm -f web_pha_${band}.gif
    if ( -e $WORK/plfit_${band}.xcm ) rm -r $WORK/plfit_${band}.xcm

    if ( -e ezyfit.result ) rm -r ezyfit.result
    if ( -e ezyfit.log ) rm -r ezyfit.log

    if ( -e ezyres_${band}.dat ) rm -r ezyres_${band}.dat
    if ( -e ezylog_${band}.dat ) rm -r ezylog_${band}.dat

end


#### copy the ezyfit_v4.tcl procedure to WORK directory ###

cp $XCO/ezyfit_v4.tcl $WORK


#### insert sequence number into scripts and run ####

# do for each band
# check the file ezyfit.result (created by the procedure ezyfit_v4.tcl) to
# see if the spectral fitting was successful

foreach band ( fullf softf hardf nhfullf nhsoftf nhhardf )

    echo "gnt_infrm (t_fits)"

    echo "gnt_infrm (t_fits) performing spectral fit for ${band}"

    cat $XCO/t_fits_pl_${band}.xcm | sed s/SEQUENCE/"$obsid"/g >! $WORK/plfit.xcm

    cp $WORK/plfit.xcm $WORK/plfit_${band}.xcm

    xspec - plfit_${band}.xcm         # run Xspec

    if ( -e ezyfit.result ) then

	set report = `tail -1 ezyfit.result`

	if ( "$#report" == "0" ) set report = "ERROR with fit"

	if ( "$report[1]" == "GOOD" ) then
	    echo "gnt_infrm (t_fits) ${band} $report"
	else
	    echo "gnt_warng (t_fits) **** ${band} $report ****"
	endif

	mv ezyfit.result ezyres_${band}.dat

    else

	echo "gnt_warng (t_fits) **** ezyfit.result was not produced for ${band} ****"

    endif


    if ( -e ezyfit.log ) then
	mv ezyfit.log ezylog_${band}.dat
    else
	echo "gnt_warng (t_fits) **** ezyfit.log was not produced for ${band} ****"
    endif

end


########################################################################################
# Check that the fit plot and log were created
########################################################################################

echo "gnt_infrm (t_fits)"
foreach band ( fullf softf hardf nhfullf nhsoftf nhhardf )
    if ( ! -e $WORK/web_pha_${band}_init.gif ) echo "gnt_warng (t_fits) **** web_pha_${band}_init.gif was not produced ****"
    if ( ! -e $WORK/web_pha_${band}.gif ) echo "gnt_warng (t_fits) **** web_pha_${band}.gif was not produced ****"
end

rm -f $WORK/xautosav.xcm     # this file is not required


########################################################################################
# Create jpeg files for the web page
########################################################################################

echo "gnt_infrm (t_fits) converting .gif web images to jpeg"

foreach band ( full soft hard )
foreach file ( web_pha_${band}f web_pha_nh${band}f web_pha_${band}f_init web_pha_nh${band}f_init )

    if ( -e ${file}.gif ) then
	echo "gnt_infrm (t_fits) creating thumbnail for file ${file}.gif"
	convert -resize 170x136 gif:${file}.gif gif:${file}_thumb.gif
    else
	echo "gnt_warng (t_fits) **** spectral image file ${file}.gif was not produced ****"
    endif

end
end


########################################################################################
# Get BACKSCAL (proportional to source region area)
########################################################################################

if ( -e backscal.dat ) rm -f backscal.dat

foreach inst ( sis0 sis1 gis2 gis3 )
    if ( -e backscal_${inst}.dat ) rm -f backscal_${inst}.dat
    set backscal = `fkeyprint ${obsid}_${inst}_src_20.pha+1 keynam=BACKSCAL | grep "BACKSCAL="`
    echo "$backscal[2]" > backscal_${inst}.dat
    echo "$backscal[2]" >> backscal.dat
end


########################################################################################
# The end of the "if detected source" statement
########################################################################################

else

    echo "gnt_warng (t_fits) **** object was not detected - no spectral fitting conducted ****"

endif


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

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