#! /bin/tcsh -f

# t_group 4.0

# Groups a pha file so that each bin has at least a set number of counts.
# the grouping in the input file is not reset, so that subsequent
# groups of bins are combined if the number of counts in the group is less
# than the required miniumum.  Technically, this is achieved by changing
# any grouping values of 1 to -1 to "combine" adjacent groups.
#
# The input spectrum is required to only have non-good (good is when QUALITY=0) channels
# at the beginning and end of the spectrum.
#
# If the last group in the spectrum has fewer than the specified minimum number of counts,
# then the channels in that group are combined with the previous group.


# Paul O'Neill

# Date in final form: 17 December 2004


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

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


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

set infile = $1           # input .pha file name
set outfile = $2          # output file name
set mincnt = $3           # minmum number of counts allowed in each group

set mincnt = `echo $mincnt | awk '{print int($1)}'`  # make sure mincnt is an integer

echo "gnt_infrm (t_group) input file name: $infile"
echo "gnt_infrm (t_group) output file name: $outfile"
echo "gnt_infrm (t_group) minimum number of counts permitted: $mincnt"


########################################################################################
# Remove any previous temporary files
########################################################################################

if ( -e t_group_indata.txt ) rm t_group_indata.txt
if ( -e t_group_outdata.txt ) rm t_group_indata.txt
if ( -e t_group_outdata.fits ) rm t_group_outdata.fits
if ( -e t_group_cols.txt ) rm t_group_cols.txt


########################################################################################
# Get values from pha file
########################################################################################


#### make sure the GROUPING column exists ####

set isgrpcol = `fstruct ${infile}+1 | grep "GROUPING"`

if ( $#isgrpcol == "0" ) then     # no GROUPING column in infile
    echo "gnt_error (t_group) file $infile has no GROUPING column"  
    fcalc infile=${infile}+1 outfile=$infile clname=GROUPING expr=1 copycol=yes copyall=yes rowrange="-" clobber=yes
    echo "gnt_infrm (t_group) GROUPING column (GROUPING=1) added to $infile" 
else
   echo "gnt_infrm (t_group) GROUPING column exists in file $infile" 
endif


#### do the fdump ####

pset fdump prhead=no prdata=yes showcol=no showunit=no showrow=yes clobber=yes
fdump infile=${infile}+1 outfile=t_group_indata.txt columns="COUNTS QUALITY GROUPING CHANNEL" rows="-"

set inrows = `awk '{print $1}' t_group_indata.txt`
set incounts = `awk '{print $2}' t_group_indata.txt`
set inqualities = `awk '{print $3}' t_group_indata.txt`
set ingroupings = `awk '{print $4}' t_group_indata.txt`
set inchannels = `awk '{print $5}' t_group_indata.txt`


#### check the lists for correct length ####

echo "gnt_infrm (t_group) number of rows in the spectrum is $#inrows"

if ( $#incounts != $#inrows ) then
    echo "gnt_fatal (t_group) the number $#incounts of counts rows is not equal to $#inrows"
    exit 1
endif

if ( $#inqualities != $#inrows ) then
    echo "gnt_fatal (t_group) the number $#inqualities of quality rows is not equal to $#inrows"
    exit 1
endif

if ( $#ingroupings != $#inrows ) then
    echo "gnt_fatal (t_group) the number $#ingroupings of grouping rows is not equal to $#inrows"
    exit 1
endif

if ( $#inchannels != $#inrows ) then
    echo "gnt_fatal (t_group) the number $#inchannels of channel rows is not equal to $#inrows"
    exit 1
endif

echo "gnt_infrm (t_group) successfully read data from $infile"
echo "gnt_infrm (t_group)"


########################################################################################
# Initialise the output lists and files 
########################################################################################

set grpqualities = ""    # used for each new group
set grpgroupings = ""

set outqualities = ""    # list of all new groups
set outgroupings = ""

@ minnumgrp = 5


########################################################################################
# Find the first GOOD and last GOOD rows and check qualtity/grouping values
########################################################################################

echo "gnt_infrm (t_group) finding first good row and last good row, and checking quality and grouping values"


#### find the first  GOOD row ####

@ row = 1

while ( ($inqualities[${row}] != 0) && ($row <= $#inrows) )
    @ row ++
end

if ( $row > $#inrows ) then
    echo "gnt_fatal (t_group) the entire spectrum is defined as BAD"
    exit 1
else
    @ firstgoodrow = $row
endif

echo "gnt_infrm (t_group) the first GOOD (QUALITY=0) row is $firstgoodrow (channel = $inchannels[${row}])"


#### find the last GOOD row ####

@ row = $#inrows

while ( ($inqualities[${row}] != 0) && ($row >= 1) )
    @ row --
end


if ( $row < 1 ) then
    echo "gnt_fatal (t_group) the entire spectrum is defined as BAD"
    exit 1
else
    @ lastgoodrow = $row
endif

echo "gnt_infrm (t_group) the last GOOD (QUALITY=0) row is $lastgoodrow (channel = $inchannels[${row}])"


#### now confirm that all "GOOD" rows are actually good ####

@ row = $firstgoodrow

while ( $row <= $lastgoodrow )
    
    if ( $inqualities[${row}] != 0 ) then
	echo "gnt_fatal (t_group) the QUALITY value is $inqualities[${row}] at row = ${row}"
	exit 1
    endif

    if ( ($ingroupings[${row}] != 1) && ($ingroupings[${row}] != -1) ) then
	echo "gnt_fatal (t_group) the GROUPING value is $ingroupings[${row}] at row = ${row}"
	exit 1
    endif

    @ row ++

end


echo "gnt_infrm (t_group) grouping and quality okay between row $firstgoodrow (channel = $inchannels[${firstgoodrow}])"
echo "gnt_infrm (t_group) and row $lastgoodrow (channel = $inchannels[${lastgoodrow}])"
echo "gnt_infrm (t_group)"
echo "gnt_infrm (t_group) proceeding with grouping"
echo "gnt_infrm (t_group)"


########################################################################################
# Step through the BAD rows at the start of the spectrum
########################################################################################

echo "gnt_infrm (t_group) stepping through BAD rows at start of spectrum"

@ row = 1

while ( $row < $firstgoodrow )
    set outqualities = ( $outqualities $inqualities[${row}] )
    set outgroupings = ( $outgroupings $ingroupings[${row}] )
    @ row ++
end


########################################################################################
# Deal with the first good row
########################################################################################

echo "gnt_infrm (t_group) dealing with first good row"

@ row = $firstgoodrow
@ grpcnt = $incounts[${row}]
set grpgroupings = 1    # start of new group
set grpqualities = 0


########################################################################################
# Step through remaining good rows
########################################################################################

echo "gnt_infrm (t_group) stepping through remaining good rows"

@ numgrp = 0     # have no new groups yet

@ row = $firstgoodrow + 1

while ( $row <= $lastgoodrow )


#### if encounter the start of a new group and have enough counts in current grp ####

    if ( $ingroupings[${row}] == "1" && $grpcnt >= $mincnt ) then
	@ numgrp ++                                        # found another group
  	set outqualities = ( $outqualities $grpqualities ) # append current group to output lists
	set outgroupings = ( $outgroupings $grpgroupings )   
	set grpqualities = ( 0 )                           # start a new list for this new group
	set grpgroupings = ( 1 )
	@ grpcnt = $incounts[${row}]                     # start counting for this new group

#### otherwise keep the grouping going ####

    else
	set grpqualities = ( $grpqualities $inqualities[${row}] )
	set grpgroupings = ( $grpgroupings -1 )
	@ grpcnt = $grpcnt + $incounts[${row}]
    endif

    @ row ++

end


########################################################################################
# Deal with the last row
########################################################################################

echo "gnt_infrm (t_group) dealing with last good row"

if ( $grpcnt >= $mincnt ) then        # final group is okay so append

    @ numgrp ++                                        # found another group
    set outqualities = ( $outqualities $grpqualities ) # append current group to output lists
    set outgroupings = ( $outgroupings $grpgroupings )
    echo "gnt_infrm (t_group)"
    echo "gnt_infrm (t_group) all rows were grouped exactly"
    echo "gnt_infrm (t_group)"

else     # this group contains spare channels - combine with last group

    foreach remrow ( $grpgroupings )
	set outqualities = ( $outqualities 0 )   # these rows are good
	set outgroupings = ( $outgroupings -1 )  # -1 to combine with previous group
    end

    echo "gnt_infrm (t_group)"
    echo "gnt_infrm (t_group) the last $#grpgroupings rows only had a total"
    echo "gnt_infrm (t_group) of $grpcnt counts, so they were grouped with the"
    echo "gnt_infrm (t_group) preceeding group"
    echo "gnt_infrm (t_group)"

endif

echo "gnt_infrm (t_group) number of good groups in output spectrum = $numgrp"
echo "gnt_infrm (t_group)"

if ( $numgrp < $minnumgrp ) then
    echo "gnt_infrm (t_group) number of good groups is less than $minnumgrp - exiting"
    exit 1
endif


########################################################################################
# Step through the BAD rows at the end of the spectrum
########################################################################################

echo "gnt_infrm (t_group) stepping through BAD rows at end of spectrum"

@ row = $lastgoodrow + 1

while ( $row <= $#inrows )
    set outqualities = ( $outqualities $inqualities[${row}] )
    set outgroupings = ( $outgroupings $ingroupings[${row}] )
    @ row ++
end


########################################################################################
# Check that output lists have the correct lengths
########################################################################################

if ( $#outqualities != $#inrows ) then
    echo "gnt_fatal (t_group) the number $#outqualities of output quality rows is not equal to $#inrows"
    exit 1
endif

if ( $#outgroupings != $#inrows ) then
    echo "gnt_fatal (t_group) the number $#outgroupings of output grouping rows is not equal to $#inrows"
    exit 1
endif



########################################################################################
# Write to file
########################################################################################


#### write to the text file ####

echo "gnt_infrm (t_group)"
echo "gnt_infrm (t_group) writing the new groupings and qualities to a text file"

touch t_group_outdata.txt    # temporary file

@ row = 1

while ( $row <= $#inrows )
    echo "$outqualities[${row}] $outgroupings[${row}]" >> t_group_outdata.txt
    @ row ++
end


#### create FITS files ####

echo "gnt_infrm (t_group) creating a FITS file with new groupings and qualities"

echo "QUALITY I" >! t_group_cols.txt
echo "GROUPING I" >> t_group_cols.txt

fcreate cdfile=t_group_cols.txt datafile=t_group_outdata.txt outfile=t_group_outdata.fits


#### create output spectrum FITS file and add columns ####

echo "gnt_infrm (t_group) creating the output FITS file"

if ( -e $outfile ) rm -f $outfile
cp $infile $outfile

fdelcol infile=${outfile}+1 colname="QUALITY" confirm=no proceed=yes
fdelcol infile=${outfile}+1 colname="GROUPING" confirm=no proceed=yes

faddcol infile=${outfile}+1 colfile=t_group_outdata.fits+1 colname="QUALITY,GROUPING"


########################################################################################
# Remove temporary files
########################################################################################

if ( -e t_group_indata.txt ) rm t_group_indata.txt
if ( -e t_group_outdata.txt ) rm t_group_outdata.txt
if ( -e t_group_outdata.fits ) rm t_group_outdata.fits
if ( -e t_group_cols.txt ) rm t_group_cols.txt

pset fdump prhead=yes prdata=yes showcol=yes showunit=yes showrow=yes clobber=yes
pset fdump outfile=STDOUT columns=""


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

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