#! /bin/tcsh -f

# t_prod 4.0

# Extracts products from the screened events files
# Uses code from do_prod (v3.2) script by K.Nandra/D.Turcan

# NOTE: do_prod_sis.xco and do_prod_gis.xco were changed from v3.2 to accept as 
# parameters the PI channels corresponding to the upper and lower boundaries
# of the soft and hard energy bands.

# Adaptation by Paul O'Neill

# Calls the scripts: do_prod_sis.xco do_prod_gis.xco

# Date in final form: 17 December 2004
 

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

echo ""
echo "gnt_infrm (t_prod) ++++++++++++++++++++++++++++++++++++++++++++++++++++++"
echo "gnt_infrm (t_prod)  t_prod 4.0  STARTED "`date`
echo "gnt_infrm (t_prod) ++++++++++++++++++++++++++++++++++++++++++++++++++++++"
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 XCO = $5             # directory with xselect/Xspec 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

echo "gnt_infrm (t_prod) running extraction procedure for $obsid"
echo " "

foreach inst ( sis0 sis1 gis2 gis3 )
    foreach ext ( .evt _src.reg _bgd.reg )
	set file = ${obsid}_${inst}$ext
	if ( ! -e $WORK/$file ) then
	    echo "gnt_fatal (t_prod) file $file does not exist."
	    exit 1
	endif
    end
end 


########################################################################################
# Fill and/or update the PI column in SIS data using sispi
########################################################################################

echo "gnt_infrm (t_prod) Running SISPI"

sispi ${obsid}_sis0.evt calfile=$CALDB/data/asca/sis/bcf/sisph2pi_290301.fits    # this is the latest file
sispi ${obsid}_sis1.evt calfile=$CALDB/data/asca/sis/bcf/sisph2pi_290301.fits    # this is the latest file


########################################################################################
# Set the hard and soft band SIS channel boundaries 
########################################################################################

# The energy-to-channel mapping shown here was determined with reference to rmf files
# created for spectra that were extracted and then saved with grouping=no.

# The energy-to-channel mapping from these response matrices is
# slightly different to that given by the equations in the ASCA ABC Guide (p 67)

if ( $datamode == "BRIGHT" ) then    # for BRIGHT mode SIS data (0.502 - 2.000 - 10.002 keV)

    set S_soft_lo = 135      # soft band, lowest PI channel  (lower = 0.502 keV)
    set S_soft2_lo = 162     # soft2 band, lowest PI channel (lower = 0.601 keV)
    set S_soft_hi = 543      # soft and soft2 band, highest PI channel (upper = 2.000 keV)

    set S_hard_lo = 544      # hard band, lowest PI channel  (lower = 2.000 keV)
    set S_hard_hi = 1707     # hard band, highest PI channel (upper = 10.002 keV)

else    # for BRIGHT2 mode SIS data (0.500 - 2.001 - 10.001 keV)

    set S_soft_lo = 134     # soft band, lowest PI channel  (lower = 0.500 keV)
    set S_soft2_lo = 162    # soft2 band, lowest PI channel (lower = 0.603 keV)
    set S_soft_hi = 543     # soft and soft2 band, highest PI channel (upper = 2.001 keV)

    set S_hard_lo = 544     # hard band, lowest PI channel  (lower = 2.001 keV)
    set S_hard_hi = 2727    # hard band, highest PI channel (upper = 10.001 keV)

endif


########################################################################################
# Set the hard and soft band GIS channel boundaries 
########################################################################################

# The energy-to-channel conversion was done with reference to the GIS2 and 
# GIS3 rmf files in CALDB
#
# GIS2: gis2v4_0.rmf
# GIS3 standard (1024 channels): gis3v4_0.rmf
# GIS3 BITFIX problem (128 channels): gis3v4_0_128ch.rmf
#
# In these rmf files the first channel number is 1.  However,
# in GIS event files the minimum channel number is 0.
# Therefore, the numbers given below are the channel numbers in the 
# rmf files, minus 1.
#
# The energy ranges for GIS2 and GIS3 data, including those data with only 128
# channels, are 0.754 - 2.057 - 10.091 keV)


#### Set channels for GIS2 ####

set G2_soft_lo = 64     # soft band, lowest PI channel  (lower energy = 0.754 keV)
set G2_soft_hi = 175    # soft band, highest PI channel (upper energy = 2.075 keV)
set G2_hard_lo = 176    # hard band, lowest PI channel  (lower energy = 2.075 keV)
set G2_hard_hi = 855    # hard band, highest PI channel (upper energy = 10.091 keV)


#### Set channels for GIS3 - may have BITFIX problem ####

if ( $bitfix == "y" ) then   # only have 128 channels

    set G3_soft_lo = 8     # soft band, lowest PI channel  (lower energy = 0.754 keV)
    set G3_soft_hi = 21    # soft band, highest PI channel (upper energy = 2.075 keV)
    set G3_hard_lo = 22    # hard band, lowest PI channel  (lower energy = 2.075 keV)
    set G3_hard_hi = 106   # hard band, highest PI channel (upper energy = 10.091 keV)

else     # set the same as GIS2 - 1024 energy channels

    set G3_soft_lo = 64     # soft band, lowest PI channel  (lower energy = 0.754 keV)
    set G3_soft_hi = 175    # soft band, highest PI channel (upper energy = 2.075 keV)
    set G3_hard_lo = 176    # hard band, lowest PI channel  (lower energy = 2.075 keV)
    set G3_hard_hi = 855    # hard band, highest PI channel (upper energy = 10.091 keV)

endif


echo "gnt_infrm (t_prod) DATAMODE = $datamode"
echo "gnt_infrm (t_prod) BITFIX = $bitfix"


########################################################################################
# Extract products using xselect
########################################################################################

foreach reg ( src bgd )

   echo "gnt_infrm (t_prod) Extracting $reg data for sis0"
   xselect @$XCO/do_prod_sis $obsid sis0 $reg 1 $S_soft_lo $S_soft_hi $S_hard_lo $S_hard_hi $S_soft2_lo

   echo "gnt_infrm (t_prod) Extracting $reg data for sis1"
   xselect @$XCO/do_prod_sis $obsid sis1 $reg 3 $S_soft_lo $S_soft_hi $S_hard_lo $S_hard_hi $S_soft2_lo

   echo "gnt_infrm (t_prod) Extracting $reg data for gis2"
   xselect @$XCO/do_prod_gis $obsid gis2 $reg $G2_soft_lo $G2_soft_hi $G2_hard_lo $G2_hard_hi 

   echo "gnt_infrm (t_prod) Extracting BITFIX=$bitfix $reg data for gis3"
   xselect @$XCO/do_prod_gis $obsid gis3 $reg $G3_soft_lo $G3_soft_hi $G3_hard_lo $G3_hard_hi 

end

########################################################################################
# Check/remove files
########################################################################################


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

foreach inst ( sis0 sis1 gis2 gis3 )

    set file = ${obsid}_${inst}_sky.img
    if ( ! -e $WORK/$file ) then
	echo "gnt_fatal (t_prod) File $file was not created."
	exit 1
    endif

    foreach ext ( .pha .evt _hard.evt _soft.evt _soft2.evt _full2.evt )
    foreach reg ( src bgd )

	set file = ${obsid}_${inst}_${reg}$ext

        if (("$inst" == "gis2" || "$inst" == "gis3") && ("$ext" == "_soft.evt" || "$ext" == "_soft2.evt" || "$ext" == "_full2.evt") ) then  
	   echo "gnt_infrm (t_prod) NOT using GIS soft, soft2, or full2 band"
	else
	    if ( ! -e $WORK/$file ) then
		echo "gnt_fatal (t_prod) File $file was not created."
		exit 1
	    endif
	endif

    end
    end

end


#### Remove unnecessary files ####

rm -f $WORK/*.tmp


########################################################################################
# Replace bad events in pixel 1,1 with a zero value
########################################################################################

fparimg 0.0 $WORK/${obsid}_sis0_sky.img 1,1 
fparimg 0.0 $WORK/${obsid}_sis1_sky.img 1,1
fparimg 0.0 $WORK/${obsid}_gis2_sky.img 1,1
fparimg 0.0 $WORK/${obsid}_gis3_sky.img 1,1


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

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