#! /bin/sh
# This is the LHEA perl script: /cvmfs/extras-fp7.egi.eu/extras/heasoft/ftools/x86_64-unknown-linux-gnu-libc2.19-0/bin/rex
# The purpose of this special block is to make this script work with
# the user's local perl, regardless of where that perl is installed.
# The variable LHEAPERL is set by the initialization script to
# point to the local perl installation.
#-------------------------------------------------------------------------------
eval '
if [ "x$LHEAPERL" = x ]; then
  echo "Please run standard LHEA initialization before attempting to run /cvmfs/extras-fp7.egi.eu/extras/heasoft/ftools/x86_64-unknown-linux-gnu-libc2.19-0/bin/rex."
  exit 3
elif [ "$LHEAPERL" = noperl ]; then
  echo "During LHEA initialization, no acceptable version of Perl was found."
  echo "Cannot execute script /cvmfs/extras-fp7.egi.eu/extras/heasoft/ftools/x86_64-unknown-linux-gnu-libc2.19-0/bin/rex."
  exit 3
elif [ `$LHEAPERL -v < /dev/null 2> /dev/null | grep -ic "perl"` -eq 0 ]; then
  echo "LHEAPERL variable does not point to a usable perl."
  exit 3
else
  # Force Perl into 32-bit mode (to match the binaries) if necessary:
  if [ "x$HD_BUILD_ARCH_32_BIT" = xyes ]; then
    if [ `$LHEAPERL -V 2> /dev/null | grep -ic "USE_64_BIT"` -ne 0 ]; then
      VERSIONER_PERL_PREFER_32_BIT=yes
      export VERSIONER_PERL_PREFER_32_BIT
    fi
  fi
  exec $LHEAPERL -x $0 ${1+"$@"}
fi
'
if(0);
# Do not delete anything above this comment from an installed LHEA script!
#-------------------------------------------------------------------------------
#!/usr1/local/bin/perl5

#  Basic script to perform extraction and background subtraction for 
#  all observations in a given proposal level directory.  

#  Uses aux directory, found in aux.tar on our ftp area.  
#  Extraction script makes file lists, GTIs, filter files, background files;
#  then extracts products for each obsid and for all together


use Getopt::Std;


getopts('d:c:l:r:p:t:x:fnsbhq');

if (defined $opt_h) {  # help flag
    print <<EOHELP1;

NAME
    rex - Makes preliminary products and then extracts data;
          i.e. creates filter files, GTIs, lists of data files, 
	  background data files and lists, lc and pha files.  


USAGE

    rex    -d "Input path to proposal level data directory"
	   -c "Channel range to extract for light curves;  0-27 by default"
	   -l "Layers to extract:  1 (default), 2, 3, or all."
           -r "Root name for output"
	   -p "Printmode:  l for light curves, s for spectra; (default both)."
           -t "Target:  give two digit target number for extraction."
           -x "To reduce HEXTE Archive mode instead, give cluster(0 1,a b)"
	   -s "Flag to process slews instead;  default is to ignore slews."
           -f "Flag to turn off filter file generation"
	   -b "Flag to turn off background generation"
           -q "Flag to quell prompting for unspecified options."
           -n "Flag to use no saahfile (for bright source models)"
	  
	   -h "Print help"

DESCRIPTION

    Rex is a script designed to run through the basic data reduction steps for
    multiple observations of a given target.  It can run with minimum input, but 
    retains the flexibility to change the defaults for most basic options.  

    The user need only give the tool the path to the FITS database proposal
    level directory;  the default procedure is then to extract both light curves 
    and spectra from layer 1 using channels 0-27 (lc's only) and PCUs 0, 1, and 2.  
    The script will create a set of subdirectories in the aux directory, one for 
    for each observation's intermediate products (file lists, GTI files, 
    background data files, logs, etc.), and another set in the working directory 
    for the resulting lightcurves and/or spectra.  For each obsid, the script will:
    
	1)  run xtefilt using the appidlist in the aux 
            directory;  -f option over-rides this;   

        2)  run maketime to create a standard list of Good Time Intervals
            using expression found in aux/expression.txt;

        4)  run pcabackest to create modeled background data for PCA
            using models found in aux/model.files and aux/pca_saa_history;
            -b option over-rides;  for HEXTE data, runs hxtback to separate
            source from offset data;

        5)  create lists of both the data and background files;

        6)  run saextrct to produce individual light curves and/or spectra;

        7)  run hxtdead to correct HEXTE spectra/light curves for deadtime;

        8)  subtract background light curve from data, and write background
            spectrum into BACKFILE keyword of data spectrum;
            
    Additionally, the script will create merged products from all the
    observations in the given directory.

    This script can be run from the command line using all default options 
    (except the data directory) by specifying the -q flag to quell the prompting.
    If the -q flag is omitted, any option not specified on the command line 
    will be prompted for, with the default taken if no value is then given. 
     
    Rex must be run from a working directory that contains the 
    "aux" directory found on the RXTE ftp site at 
      ftp://legacy.gsfc.nasa.gov/xte/software/rex/aux.tar
    In this tar file are the miscellaneous files that will be used by the script, 
    namely the goodtime expression used in maketime, the column lists to be 
    extracted, and the background model files.  The text files can be edited, but
    the names of the text files must remain the same in order for the script 
    to find them.  

    While this script is designed for ease of use rather than flexibility, 
    the existence of several optional inputs and the ability to  
    edit the tools in the aux directory give the user some control.  To use 
    use different model files, edit the model.files text list to point
    to the new files;  to change the number of PCUs extracted, add to the column
    list files, etc. Time consuming steps which may have already been performed 
    such as xtefilt and pcabackest can be bypassed using the flags -f and -b 
    respectively.  Changes to the channel boundaries and layers can be 
    made using the -c and -l options, respectively.  To limit extraction 
    to only lightcurves or only spectra, us the -p option.  Changes to the 
    criteria for determining Good Time Intervals (GTIs) can be made by 
    editing the expression in aux/expression.txt.  By default, the 
    script will only extract data from PCUs 0, 1, and 2, as PCUs 3 and 
    4 usually go on and off in the course of a long series of 
    observations.  To add PCUs 3 and 4 to the extraction, the user can 
    edit the column files in the aux directory to add the additional 
    columns.  (The selection expression should then be modified to include
    a test to make sure these two PCUs are on.)

    In a script this long, the possibilities for errors are numerous.  Most should
    be trapped by the script, which will give a warning, write any error messages
    from the FTOOLS into a log in the aux directory, and either continue if possible
    or quit.  When a warning appears, check the log file for the error message from
    the FTOOL;  if this message does not enlighten, send the log file on to 
    xtehelp\@athena.gsfc.nasa.gov along with the inputs you gave to Rex 
    originally and any information about what you may have changed in the 
    aux directory.
    


PARAMETERS (SWITCHES/FLAGS)

          -d [path]  - Input the path to the proposal level of your FITS 
                       database, e.g. /local/data/P30402. 

          -c [l-u]   - Give the desired channel range for the resulting light 
                       curves, e.g. 5-30.  If nothing is specified, 0-27 
                       will be extracted.

          -l [layer] - Specify the layer to be extracted:  1, 2, 3, or all 
                       are the possible entries.  Anything else will be 
                       ignored and the default used.  If nothing is 
                       specified, layer 1 only will be extracted.  

          -r [root]  - Root name for output.  Files will be named with this 
                       root, if specified, followed by layerN, chX-Y (lc's 
                       only), and bkg for background files.  

          -p [mode]  - Specify printmode, l for light curves only, s for 
                       spectra only.  Default is both.  

          -t [targ]  - Give two digit target number, e.g. "03".  Default 
                       is all found.

          -x [cl]    - To reduce HEXTE Archive mode data instead of PCA Standard2,
                       give the cluster you wish to analyze.  This option works
                       with all options except the layer option.

          -s         - Flag to include processing of slew data.  By default, 
                       the  script will ignore slew observations, i.e. those 
                       whose ObsIds end in A, Z, C, or D.

          -f         - Flag to skip filter file regeneration.  If the filter 
                       files have already been made and include all the 
                       necessary columns in aux/expression.txt and those 
                       needed by pcabackest (PCUn_ON, BKGD_THETA, BKGD_PHI), 
                       this time consuming step can  be avoided by specifying 
                       this flag.  The filter files must exist in the FITS 
                       database, in the stdprod directories for each obsid, 
                       and they must be refered to by the FIST index file.  
                       (If made by xtefilt with the -s option, all should be 
                       well.)

          -b         - Flag to bypass the making of the background files.  
                       Use this option only if the script has already been 
                       run successfully  and you simply want to re-extract 
                       light curves and spectra.  

          -n         - Flag to use no SAA history file in pcabackest;  this 
                       is only for use with other models than the Faint L7 
                       models.  It is only necessary when there is no current
                       SAA history for the data.  
  
          -q         - Flag to quell prompting for options not specified on the
                       command line.  Anything not specified will be set to the 
                       default.

EXAMPLES

    1)  rex

        Prompts for data, target, printmode, channel, layer, and root name.  

    2)  rex -d /local/data/P30402 -q

        Runs using all defaults.

    3)  rex -bfd /local/data/P30402 -p l -l 3 -c 30-50 -q

        Runs without generating new filter files, using existing background
        files, printing only light curves, and extracting only layer 3 and
        channels 30-50.  

SEE ALSO

        SAEXTRCT, PCABACKEST, XTEFILT, MAKETIME, LCMATH, HXTBACK, HXTDEAD, REBINLC, MATHPHA 

EOHELP1


exit;
    }

$| = 1;   ## "print buffer flag" to print everything immediately, rather 
          ##  than store output until it's ready

print "\n************************************************\n";
print "    Rex script for RXTE Data Reduction v.0.3\n";
print "************************************************\n\n";

die "Directory aux not found;  please read fhelp and then get and untar \nthe auxiliary files Rex needs from: \n\n ftp://legacy.gsfc.nasa.gov/xte/software/rex/aux.tar\nRex quitting.\n" unless (-d "aux");


#############################################################################
## PRELIMINARIES:
## get command line options or prompt
#############################################################################

# get data dir
if (defined $opt_d) {
    $data = $opt_d;
}
else {
    print "\nPath to proposal level directory containing FMI: ";
    chomp($data=<STDIN>);
    $data =~ s/\s+//g;
}

# get target (OK for $target not defined)
if (defined $opt_t) {
    $target = $opt_t;
}
elsif (!defined($opt_q)) {
    print "Target [all]: ";
    chomp($target=<STDIN>);
    $target =~ s/\s+//g;
    $target = "" unless ($target =~ /\d\d/);
}

#get printmode (OK if $opt_p undefined)
if (!defined($opt_p) && !defined($opt_q)) {
    print "Printmode, l (ligh curve), s (spectra) [both]: ";
    chomp($opt_p = <STDIN>);
    $opt_p =~ s/\s+//g;
    $opt_p = "both" unless ($opt_p);
}

#get channels
if (defined $opt_c) {
    $channels = $opt_c;
}
elsif (!defined($opt_q)) {
    print "Channels for light curve extraction [0-27]: ";
    chomp($channels = <STDIN>);
    $channels =~ s/\s+//g;
    $channels = "0-27" unless ($channels =~ /\d+-\d+/);
}
else {
    $channels = "0-27";
}

# get layer
if ((defined $opt_l) && !defined($opt_x)) {
  # if PCA and layer is given correctly, use it
    if ($opt_l =~ /^[123]$/) {
	$colfile = "layer${opt_l}.cols";
        $layer = $opt_l;
    }
    elsif ($opt_l eq "all") {
	$colfile = "allayers.cols";
        $layer = $opt_l;
    }
    else {
	$colfile = "layer1.cols";
        $layer = "1";
    }
}
elsif (!defined($opt_q) && !defined($opt_x)) {
  # if PCA and not quelled, then prompt
    print "Layers to be extracted (1,2,3, or all) [1]: ";
    chomp($layer = <STDIN>);
    $layer =~ s/\s+//g;
    if (!$layer) {
	$colfile = "layer1.cols";
	$layer = "1";
    }
    elsif ($layer =~ /^[123]$/) {
	$colfile = "layer${layer}.cols";
    }
    elsif ($layer eq "all") {
	$colfile = "allayers.cols";
    }
    else {
	$colfile = "layer1.cols";
        $layer = "1";
    }

}
elsif (!defined($opt_x)) { 
  # if PCA and quelled, use layer 1
    $layer = "1";
    $colfile = "layer1.cols";
}
else {
  # otherwise use basic HEXTE column list
  $colfile = "hexte.cols";
}

# get root name
if (defined($opt_r)) {
    $root = "${opt_r}";
}
elsif (!defined($opt_q)) {
    print "Root name for products []: ";
    chomp($root = <STDIN>);
    $root =~ s/\s+//g;
}

##
## print out inputs/options
##

if (-e "${data}/FMI") {
    print "\nDATA:\t\t$data\n";
}
else {
    die "\nFMI not found in ${data}.\nPlease give path to proposal level directory containing FMI.\n";
}
print "TARGET:\t\t$target\n" if ($target);
print "TARGET:\t\tall\n" unless ($target);

if ($opt_p eq "l") {
    print "PRINTMODE:\tlight curves\n";
}
elsif ($opt_p eq "s") {
    print "PRINTMODE:\tspectra\n";
}
else {
    print "PRINTMODE:\tboth\n";
}

print "CHANNELS:\t$channels\n";

print "LAYER:\t\t$layer\n" unless (defined($opt_x));

print "ROOT NAME:\t$root\n" if ($root);

$root .= "_" if ($root);

#this is the change to make Rex accept a and b's (resets them to 0 and 1)
if (defined($opt_x)) {
  if ($opt_x =~ /a|A/) {
    $opt_x = 0;
  }
  if ($opt_x =~ /b|B/) {
    $opt_x = 1;
  }
}

# set output filenames using root

if (defined($opt_x)) {
  $pharoot = "$root"."hxt${opt_x}";
  $lcroot = "$root"."hxt${opt_x}_ch${channels}";
}
else {
  $pharoot = "$root"."layer${layer}";
  $lcroot = "$root"."layer${layer}_ch${channels}";
}

$expr = `cat aux/expression.txt`;
chomp $expr;
print "EXPRESSION:\t$expr\n\n";

if (defined $opt_s) {
    print "Using slew observations only.\n";
}

print "Making filter files for all obsids in $data \n" unless (defined($opt_f));
print "Skipping filter file generation. \n" if (defined($opt_f));

print "Skipping background generation. \n" if (defined($opt_b) && !defined($opt_x));
print "Skipping hxtback step. \n" if (defined($opt_b) && defined($opt_x));

print "Analayzing HEXTE Archive mode data for cluster $opt_x\n" if (defined($opt_x));

print "\n";

#############################################################################
## SETUP:
## fdump the FMI file and open full lists of PCA Std2 or HEXTE archive data, 
## filter files, and bkgd; 
## initialize everything
#############################################################################

## define what you're doing, PCA or HEXTE

if (defined($opt_x)) { ## define HEXTE 
  $indat = "hxtsrc";
  $bakdat = "hxtbkg";
  $sys = "hexte";
}
else {  ## define PCA
  $indat = "std2";
  $bakdat = "bkg";
  $sys = "pca";
}

## dump appropriate FMI columns
$fmistring = &fdmp("${data}/FMI","obsid PCA_Index_File Std_Prod_Index_File StopDate") unless (defined($opt_x));
$fmistring = &fdmp("${data}/FMI","obsid HEXTE_Index_File Std_Prod_Index_File") if (defined($opt_x));

if ($fmistring) {
    @fmiline = split('\n',$fmistring);
}
else {
    die "\nProblem with ${data}/FMI:  $!\nRex quitting.\n";
}

$path = `pwd`;  ## storing path for prepending to ${bakdat}d file lists
chomp $path;

$saexpar = " gtiorfile=APPLY accumulate=ONE timecol=TIME binsz=16 lcmode=RATE spmode=SUM timemin=INDEF timemax=INDEF timeint=INDEF chmin=INDEF chmax=INDEF  chbin=INDEF bailout=yes clobber=yes ";


################################################################################
##  MAIN LOOP through obsids
################################################################################

foreach $line (@fmiline){  #  for each row in FMI, get FI and FIST
  if ($line =~ /\w/){  ## for each (non-blank) line, i.e. obsid in FMI

    @entries = split(' ',$line);
    $obsid = $entries[0];

    ##
    ## testing
    ##

    ## skip this one if the obsid is a slew and -s not flagged
    next if (!defined($opt_s) && ($obsid =~ /[AZCD]$/));

    ## skip this one if the obsids is not a slew and -s is flagged
    next if (defined($opt_s) && !($obsid =~ /[AZCD]$/));

    ## skip this one if it's not the specified target
    next if ( $target && !($obsid =~ /^\d\d\d\d\d-$target/) );
    
    if (!(-d "${data}/$obsid")) {  ## run only if obsid exists
      print "Obisd $obsid not present in $data;  skipping.\n";
      next;
    }

    $fi = $entries[1];
    $fist = $entries[2];

    ##
    ##  Check date of Observation versus SAA history validity
    ##

    if (!(defined($opt_x)) && !(defined($opt_b))) {  

      $obsend = $entries[3];
      $obsend =~ s/(\d{2})\/(\d{2})\/(\d{2})/$3$2$1/g;
      # MJT 19Jan2000 -- pivot on 1990 for Y2K!
      if ($3 > 90) {
	$obsend = $obsend + 19000000
      } else {
	$obsend = $obsend + 20000000
      }
      
      $saah = `fkeyprint aux/pca_saa_history CVED0001 | grep =`;
      $saah =~ s/^.*'(\d{2})(\d{2})-(\d{2})-(\d{2})'.*$/$1$2$3$4/g;
      chomp($saah);
      
      if ($saah < $obsend) {
	
	print "SAA history file does not cover obsid $obsid;  skipping.\n";
	next;
      }
    }

    ##
    ##  Make Std2 file list and run xtefilt for each obsid.  
    ##

    if (-e "${data}/${fi}") {
      system("mkdir aux/${obsid}") unless (-e "aux/${obsid}");
    }
    else {
      print "Cannot find subsytem index file for $obsid;  skipping\n";
      next;
    }
    # dump appropriate columns from subsystem index file
    $fidump = &fdmp("${data}/${fi}","EA6Data") unless (defined($opt_x));
    $fidump = &fdmp("${data}/${fi}","Cl${opt_x}ArchData") if (defined($opt_x));
    @fijunk = split('\s*\n',$fidump);

    ##
    ## make a hash using the file name as key, which will overwrite  
    ## redundant terms;  then print them in order in the ${indat}.xdf file
    ## in aux/${obsid} using the path ${data}/${obsid}/${sys} 
    ##
    foreach (keys(%indathash)) {  ## remove entries in ${indat} hash from
      delete $indathash{$_};     ## previous obsid runs
    }

  foreach (@fijunk) {  
    @temp = split('/');
    $indathash{$temp[1]}++ if $_;
  }
  if (!(defined($opt_x))) { ## make data file list for PCA
    unlink "aux/${obsid}/${indat}.xdf";
    open(DATA,">>aux/${obsid}/${indat}.xdf") or die "\nCannot open aux/${obsid}/${indat}.xdf to write:  $! \nRex quitting.\n";

    foreach (sort(keys(%indathash))) {
      print DATA "${data}/${obsid}/pca/$_\n";
    }
    close(DATA);
  }
    # for HEXTE, wait until after hxtback to write out any data file names
    # instead, get the 16s Hk file
  else {
    $hkjunk = &fdmp("${data}/${fi}","Cl${opt_x}Hk_16");
    @junk = split(' ',$hkjunk);
    $hk16 = $junk[0]; # first non-empty row entry
  }

    ##
    ## run xtefilt for the obsid, placing in stdprod directory
    ## then dump index file to get current filter
    ## (skips xtefilt if -f flagged for filters already made)
    ##

    if (!defined($opt_f)) {  ## unless flagged to skip xtefilt

      $args = "xtefilt -s -c -a aux/appidlist -o $obsid -p $data -t 16.0";
      unlink "aux/${obsid}/filt.log";
      open(LOG,">aux/${obsid}/filt.log");
      print "Running xtefilt for $obsid \n";
      @log = &call($args,"warn","Problem running:\n $args","aux/${obsid}/filt.log");

      $filtstat = shift(@log);
      print LOG @log;
      close(LOG);
      ## check status of filer file run;  skip obsid if error
      if ($filtstat) {
	print "Skipping obsid $obsid\n";
	next;
      }

    } # end of filter file 

    ## dump FIST file to get filter file name for maketime call 
    ## and merge list

    if (-e "${data}/${fist}") {
      $fistdump = &fdmp("${data}/${fist}","Filter_File");      
      @fistjunk = split(' ',$fistdump);
      $fistjunk[0] =~ s/stdprod\///;
      $filter = $fistjunk[0];
    }
    else {
      print "Cannot find Stdprod Index File for $obsid;  skipping\n";
      next;
    }

    ##
    ## loop through ${indat} files, running pcabackest or hxtback on each, and entering
    ## resulting background file in aux/obsid/${bakdat}.xdf 
    ## only if $opt_b not defined;  will skip if flagged
    ##


    if (!defined($opt_x) && !(defined($opt_b))) { ## PCA background generation

      print "Running pcabackest on all standard2 files in $obsid\n";
      unlink "aux/${obsid}/${bakdat}.xdf";
      open(BKG,">>aux/${obsid}/${bakdat}.xdf") or die "\nCannot open ${bakdat}.xdf to write:  $!\nRex quitting.\n";
      unlink "aux/${obsid}/${bakdat}.log";
      open(LOG,">>aux/${obsid}/${bakdat}.log") or die "\nCannot open ${bakdat}.log to write:  $!\nRex quitting.\n";

      foreach (sort(keys(%indathash))) {
	# run pcabackest for each standard2 data file
	$std2file = $_;
	$bkgfile = "aux/${obsid}/${std2file}_bkg"; 
	$args = "pcabackest infile=${data}/${obsid}/pca/${std2file} outfile=${bkgfile} ";
	$args .= " filterfile=${data}/${obsid}/stdprod/${filter} ";
	$args .= " modelfile=\@aux/model.files ";
	$args .= " modeltype=both interval=16";
	$args .= " saahfile=aux/pca_saa_history" unless (defined($opt_n));
	$args .= " saahfile=none" if (defined($opt_n));
	$args .= " propane=no layers=yes gaincorr=no ";
	$args .= " gcorrfile=- fullspec=no  syserr=no clobber=yes timeslop=128 ";
	
	@log = &call($args,"warn","Problem running:\n $args","aux/${obsid}/bkg.log");

	##
	## check status of pcabackest before writing the resulting file into the lists
	##
	print LOG "@log \n";
	if (!shift(@log)) {  ## check error status;  if 0, add to lists
	  print LOG "@log \n";
	  print BKG "${path}/${bkgfile}\n";
	} 

      }  ##  end of loop through std2.xdf files
      close(BKG);
      close(LOG);
    }  ##  end of PCA bkgd loop if $opt_b not defined

    elsif (!(defined($opt_b))) {  
    ## run hxtback on all input HEXTE FILES

      print "Running hxtback on all HEXTE archive mode files in $obsid\n";
      unlink "aux/${obsid}/${bakdat}.xdf";
      open(BKG,">>aux/${obsid}/${bakdat}.xdf") or die "\nCannot open ${bakdat}.xdf to write:  $!\nRex quitting.\n";
      unlink "aux/${obsid}/${bakdat}.log";
      open(LOG,">>aux/${obsid}/${bakdat}.log") or die "\nCannot open ${bakdat}.log to write:  $!\nRex quitting.\n";
      
      unlink "aux/${obsid}/${indat}.xdf";
      open(DATA,">>aux/${obsid}/${indat}.xdf") or die "\nCannot open aux/${obsid}/${indat}.xdf to write:  $! \nRex quitting.\n";

      foreach (sort(keys(%indathash))) {
	# for each input data file, run hxtback and enter src and bkg files in lists
	$file = $_;
	if (!(-e "${path}/aux/${obsid}/${file}_src")) {
	  # check that it hasn't already been run (it doesn't fail the way pcabackest might;  can't see why you'd want to overwrite an existing version)
	  $args = "hxtback -b -i ${data}/${obsid}/hexte/$file";
	  @log = &call($args,"warn","Problem running:\n $args","aux/${obsid}/hxtbkg.log");
	  print LOG "@log\n";
	  $stat = shift(@log);
	  if (!($stat)) {
	    # since hxtback writes output here, move into aux/obsid
	    system("mv ${file}_src ${path}/aux/${obsid}/");
	    system("mv ${file}_bkg ${path}/aux/${obsid}/");
	  }
	} # end of if not already run
	# lists recompiled anyway
	print DATA "${path}/aux/${obsid}/${file}_src\n";
	print BKG "${path}/aux/${obsid}/${file}_bkg\n";
      } ## end of hxtback run
      close(DATA);
      close(BKG);
      close(LOG);
    } # end of hxtback

    ##
    ## run maketime for individual obsids, each with log aux/obsid/maketime.log
    ## 

    $args = "maketime infile=${data}/${obsid}/stdprod/${filter} outfile=aux/${obsid}/basic.gti expr=\"${expr}\" name=NAME value=VALUE time=TIME compact=no clobber=yes ";
    unlink "aux/${obsid}/basic.gti";

    @log = &call($args,"warn","Problem running:\n$args\n\tSkipping extraction of obsid $obsid.\n","aux/${obsid}/maketime.log");
    $gtistat = shift(@log);
    ## write log only if there is an error
    if ($gtistat) {
      print "Error running maketime on $obsid;  skipping\n";
      unlink "aux/${obsid}/maketime.log";
      open(LOG,">aux/${obsid}/maketime.log") or die "\nCannot open aux/${obsid}/maketime.log to write:  $!\nRex quitting.\n";
      print LOG @log;
      close(LOG);
      next;
    }

    ## test for rows before extracting;  if not, skip to next obsid with warning
    $gti = &fdmp("aux/${obsid}/basic.gti","-");
    if (!($gti =~ /\w/)) {
      print "Observation $obsid has no good time using this expression; \n\t skipping extraction.\n";
      next;
    }
    ## test for good bgfiles
    $bgnum = (-s "aux/${obsid}/${bakdat}.xdf");  
    print "No background files produced sucessfully for $obsid; \n\t skipping extraction of background.\n" unless ($bgnum);

    # create output obsid directory
    system("mkdir $obsid") unless (-e "$obsid");

    print "Running saextrct for $obsid\n" unless (defined($opt_x));
    print "Running saextrct, hxtdead, and merging for $obsid\n" if (defined($opt_x));
    unlink "aux/${obsid}/extract.log";
    open(LOG,">>aux/${obsid}/extract.log") or die "\nCannot open aux/${obsid}/extract.log to write:  $!\nRex quitting.\n";
    
    if (!defined($opt_x)) {

    ############################################################################
    ## PCA extraction
    ## run saextrct for individual obsids, each with log aux/obsid/extract.log
    ############################################################################

      if (!($opt_p eq "l")) {  ## unless flagged for lc only

	$args = "saextrct infile=\@aux/${obsid}/${indat}.xdf gtiandfile=aux/${obsid}/basic.gti outroot=${obsid}/${pharoot} columns=\@aux/${colfile}  printmode=SPECTRUM chint=INDEF ";
	$args .= $saexpar;

	@log = &call($args,"warn","Problem running:\n $args","aux/${obsid}/extract.log");
	$phstat1 = shift(@log);
	$allphas .= "${obsid}/${pharoot}.pha " unless ($phstat1);
	print LOG @log;

	if ($bgnum) {
	  $args = "saextrct infile=\@aux/${obsid}/${bakdat}.xdf gtiandfile=aux/${obsid}/basic.gti outroot=${obsid}/${pharoot}_bkg columns=\@aux/${colfile} printmode=SPECTRUM chint=INDEF";
	  $args .= $saexpar;
	  
	  @log = &call($args,"warn","Problem running:\n $args","aux/${obsid}/extract.log");
	  $phstat2 = shift(@log);
	  print LOG @log;
	    
	  $args = "fparkey value=${obsid}/${pharoot}_bkg.pha fitsfile=${obsid}/${pharoot}.pha keyword=BACKFILE";

	  ## added check that not only are there no errors (phstat) but also
	  ## that the products exhist;  saextrct can die very quietly
	  system($args) unless (($phstat1 + $phstat2) || !((-e "${obsid}/${pharoot}.pha") && (-e "${obsid}/${pharoot}_bkg.pha")));
	  $allphasbkg .= "${obsid}/${pharoot}_bkg.pha " unless ($phstat2);

	} ## if bgnum good

      }  ## spectra 

      if (!($opt_p eq "s")) { ## unless flagged for spectra only
	
	$args = "saextrct infile=\@aux/${obsid}/${indat}.xdf gtiandfile=aux/${obsid}/basic.gti outroot=${obsid}/${lcroot} columns=\@aux/${colfile} printmode=LIGHTCURVE chint=${channels} ";
	$args .= $saexpar;
	@log = &call($args,"warn","Problem running:\n $args","aux/${obsid}/extract.log");
	$lcstat1 = shift(@log);
	print LOG @log;
	$alllcs .= "${obsid}/${lcroot}.lc " unless ($lcstat1);
	    
	if ($bgnum) {
	  $args = "saextrct infile=\@aux/${obsid}/${bakdat}.xdf  gtiandfile=aux/${obsid}/basic.gti outroot=${obsid}/${lcroot}_bkg  columns=\@aux/${colfile} printmode=LIGHTCURVE chint=${channels} ";
	  $args .= $saexpar;

	  @log = &call($args,"warn","Problem running:\n $args","aux/${obsid}/extract.log");
	  $lcstat2 = shift(@log);
	  print LOG @log;
	  close(LOG);
	  $alllcsbkg .= "${obsid}/${lcroot}_bkg.lc " unless ($lcstat2);
	  
	  ## subtract background lc from data using lcmath
		
	  $args = "lcmath infile=${obsid}/${lcroot}.lc bgfile=${obsid}/${lcroot}_bkg.lc outfile=${obsid}/${lcroot}_net.lc multi=1 multb=1 addsubr=no ";
	  
	  ## because lcmath won't clobber
	  unlink "${obsid}/${lcroot}_net.lc";
	  ## added check that not only are there no errors (lcstat) but 
	  ## also that the products exhist;  saextrct can die quietly
	  @log = `$args` unless (($lcstat1 + $lcstat2) || !((-e "${obsid}/${lcroot}.lc") && (-e "${obsid}/${lcroot}_bkg.lc")));
	  
	}  ## if bgnum good
	
      }  ## unless flagged for spectra

    } ## end of PCA extraction

    ############################################################################
    ## HEXTE extraction
    ## run saextrct for individual obsids, each with log aux/obsid/extract.log
    ############################################################################
    else {
      unlink("aux/${obsid}/hxtdead.log");
      open(DLOG,">>aux/${obsid}/hxtdead.log") or die "Cannot open aux/${obsid}/hxtdead.log to write:$!\n";

      # same code for source and background
      foreach $dat ("src","bkg") {
	open(DATA,"aux/${obsid}/hxt${dat}.xdf") or die "Cannot open aux/${obsid}/hxt${dat}.xdf to read:  $!\n";
	##
	## looping through data files one at a time because hxtdead only handles one 
	##  at a time;  extract, deadtime correct, merge results
	##
	$n = 0;
	$lclist = "";
	$phalist = "";
	while(<DATA>) {
	  $file = $_;
	  chomp $file;
	  if (!($opt_p eq "l")) {
	    $args = "saextrct infile=$file gtiandfile=aux/${obsid}/basic.gti outroot=${dat}_${n} columns=\@aux/${colfile} printmode=SPECTRUM chint=INDEF ";
	    $args .= $saexpar;
	    @log =  &call($args,"warn","Problem running:\n $args","aux/${obsid}/extract.log");
	    $phstat1 = shift(@log);

	    print LOG @log;
	    if (!$phstat1 && (-e "${dat}_${n}.pha")) {
	      $args = "hxtdead calvalf=caldb engvalf=${data}/${obsid}/$hk16 eventarcf=$file phalcf=${dat}_${n}.pha";
	      @log = &call($args,"warn","Problem running:\n $args\n","aux/${obsid}/hxtdead.log");
	      print DLOG @log;
	      $phalist .= "${dat}_${n}.pha ";
	    } # test
	  } # end of source spectra
	  if (!($opt_p eq "s")) {
	    $args = "saextrct infile=$file gtiandfile=aux/${obsid}/basic.gti outroot=${dat}_${n} columns=\@aux/${colfile} printmode=LIGHTCURVE chint=${channels} ";
	    $args .= $saexpar;
	    @log =  &call($args,"warn","Problem running:\n $args","aux/${obsid}/extract.log");
	    $lcstat1 = shift(@log);
	    print LOG @log;
	    if (!$lcstat1 && (-e "${dat}_${n}.lc")) {
	      $args = "hxtdead calvalf=caldb engvalf=${data}/${obsid}/$hk16 eventarcf=$file phalcf=${dat}_${n}.lc";
	      @log = &call($args,"warn","Problem running: \n$args\n","aux/${obsid}/hxtdead.log");
	      print DLOG @log;
	      $lclist .= "${dat}_${n}.lc ";
	    } # test
	  } # end of source lightcurves
	    $n++;
	} # while <DATA>
	
	# fmerge and lcmath lightcurves, mathpha spectra etc.

	if (!($opt_p eq "l") && (scalar(split(' ',$phalist)) > 0) ) {

	  if (scalar(split(' ',$phalist)) > 1) {
	    
	    # merge spectra
	    unlink("aux/${obsid}/mathpha.log");
	    &mergepha("$phalist","${obsid}/${pharoot}_bkg.pha","aux/${obsid}/mathpha.log") if ($dat eq "bkg");
	    &mergepha("$phalist","${obsid}/${pharoot}.pha","aux/${obsid}/mathpha.log") if ($dat eq "src");
	  }
	  elsif (scalar(split(' ',$phalist)) == 1) {
	    # or rename one spectrum 
	    system("mv $phalist ${obsid}/${pharoot}.pha") if ($dat eq "src");
	    system("mv $phalist ${obsid}/${pharoot}_bkg.pha") if ($dat eq "bkg");
	  }
	  # add this obsids spectra to master list for merging at the end
	  $allphas .= "${obsid}/${pharoot}.pha " if ($dat eq "src");
	  $allphasbkg .= "${obsid}/${pharoot}_bkg.pha " if ($dat eq "bkg");
	  # remove individual spectra
	  unlink split(' ',$phalist);
	}
	if (!($opt_p eq "s") && (scalar(split(' ',$lclist)) > 0 ) ) { 

	  if (scalar(split(' ',$lclist)) > 1 ) {
	    # merge light curves
	    @log = &mergelcs("$lclist","${obsid}/${lcroot}_bkg.lc","aux/${obsid}/merge.log") if ($dat eq "bkg");
	    @log = &mergelcs("$lclist","${obsid}/${lcroot}.lc","aux/${obsid}/merge.log") if ($dat eq "src");

	    if (shift(@log)) {
	      open(LOG,">aux/${obsid}/merge.log") or die "Cannot open aux/${obsid}/merge.log to write:  $!\n";
	      print LOG "@log";
	      close(LOG);
	      next;
	    }
	  }
	  elsif (scalar(split(' ',$lclist)) == 1 ) {
	    # or rename one light curve 
	    system("mv $lclist ${obsid}/${lcroot}.lc") if ($dat eq "src");
	    system("mv $lclist ${obsid}/${lcroot}_bkg.lc") if ($dat eq "bkg");
	  } # end of if only one

	  $alllcs .= "${obsid}/${lcroot}.lc " if ($dat eq "src");
	  $alllcsbkg .= "${obsid}/${lcroot}_${dat}.lc " if ($dat eq "bkg");	  
	  # remove individual light curves
	  unlink split(' ',$lclist);

	} # end of merging light curves

      } # end of foreach src and bkg

      # rebin background lightcurve to 64 seconds so it can be subtracted from source
      if (!($opt_p eq "s") && (-e "${obsid}/${lcroot}_bkg.lc")) {
	$args = "rebinlc infile=${obsid}/${lcroot}_bkg.lc multiple=4 outfile=\!${obsid}/${lcroot}_bkg.lc";
	@log = `$args`;
	# ... and subtract
	$args = "lcmath infile=${obsid}/${lcroot}.lc bgfile=${obsid}/${lcroot}_bkg.lc outfile=${obsid}/${lcroot}_net.lc multi=1 multb=1 addsubr=no";
	unlink("${obsid}/${lcroot}_net.lc");
	@log = `$args`;
      }
      if (!($opt_p eq "l") && (-e "${obsid}/$pharoot}_bkg")) {
	$args = "fparkey value=${obsid}/${pharoot}_bkg.pha fitsfile=${obsid}/${pharoot}.pha keyword=BACKFILE";
	system("$args") unless (($phstat1 + $phstat2) || !((-e "${obsid}/${pharoot}.pha") && (-e "${obsid}/${pharoot}_bkg.pha")));
      }
    } # else end of hexte
    
  }  ## if the line isn't blank
}  ## foreach line in FMI dump, i.e for each obsid


#############################################################################
##  End of loop through obsids.  Now, 
##  Make merged products for all obsids together
#############################################################################

if (scalar(@fmiline) == 1) {  ## only one obsid ( @fmiline has 3 returns before 1st entry)
    print "Rex finished.\n";
    exit;
}

## first, test for #files < 100 and # obsids > 1

$wc = scalar(split(' ',$allphas)) unless ($opt_p eq "l");
$wc = scalar(split(' ',$alllcs)) if ($opt_p eq "l");
if ($wc == 0) {
    print "No obsids found.\n";
    print "Rex finished.\n";
    exit;
}

print "Merging all obsids\n"; # old 100 file limit is avoided by using sumpha

print "No good background files created;  skipping background for all obsids.\n" unless (scalar(split(' ',"$alllcsbkg")) or (scalar(split(' ',"$allphasbkg"))) );

##
## merge spectra
##

if (!($opt_p eq "l")) {  
## unless flagged for lc's only

  # merge
  unlink("aux/mathpha.log");
  &mergepha("$allphas","${pharoot}.pha","aux/mathpha.log");

  ## same for background spectra
  if (scalar(split(' ',"$allphasbkg"))) { 

      &mergepha("$allphasbkg","${pharoot}_bkg.pha","aux/mathpha.log");
      
      # insert background spectrum into BACKFILE keyword of source spectrum
      if ((-e "${pharoot}_bkg.pha") and (-e "${pharoot}.pha")){
	$args = "fparkey value=${pharoot}_bkg.pha fitsfile=${pharoot}.pha keyword=BACKFILE";
	system($args);
	$args = "fchecksum infile=${pharoot}.pha update=yes datasum=yes";
	system($args);
      }
      
    }  ##  background files 
}  ## end of spectra

##
## merge light curves
##

if (!($opt_p eq "s")) {  ## unless flagged for pha's only
  # fmerge all the light curves
  @log = &mergelcs("$alllcs","$lcroot.lc","aux/merge.log");
  if (shift(@log)) {
    open(LOG,">>aux/merge.log") or die "Cannot open aux/merge.log to write:  $!\n";
    print LOG @log;
    close LOG;
  }
  
  if (scalar(split(' ',"$alllcsbkg"))) {  ## background
    @log = &mergelcs("${alllcsbkg}","${lcroot}_bkg.lc","aux/merge.log");
    
    if (shift(@log)) {
      open(LOG,">>aux/merge.log") or die "Cannot open aux/merge.log to write:  $!\n";
     print LOG @log;
      close LOG;
    }
    else {
         
      $args = "lcmath infile=${lcroot}.lc bgfile=${lcroot}_bkg.lc outfile=${lcroot}_net.lc multi=1 multb=1 addsubr=no";
      unlink("${lcroot}_net.lc");
      @log = `$args` if ((-e "${lcroot}.lc") && (-e "${lcroot}_bkg.lc"));
    
    }
  }  ## background
}  ## end of light curves

close(LOG);

print "Rex finished.\n";

#################################################################################
#  END OF MAIN 
#################################################################################


sub call {

## based on runcom, but keeps the output from the ftool even in the case
## of an error;  give command, error routine, error string, error log;
## returns array:  1st value is 0 if good or $? or 1 if system or ftool error; 
## then command; then error messages sent by ftool or system to 
## STDERR; then any other output from tool;  retains suppress option.

    local($command) = $_[0];
    local($err_routine) = $_[1];
    local($err_string) = $_[2];
    local($err_log) = $_[3];
    local($ERRFIL,@result,$CURHANDLE,$errormsg);

    if($command =~ s/^-// ) {local($suppress) = 1;}
    
    $ERRFIL = "/tmp/error$$";
    unlink $ERRFIL;
    $command .= " 2> $ERRFIL";
    @result = `$command`;

    if ( ! open(ERRFIL) ) {  ## if no ERRFIL created, return with $? status
        unlink $ERRFIL;
    }
    else {
      while ( ($_ = <ERRFIL>) =~ /(ld\.so|^\s*\n$)/ ) { ;} 
      if ( ! $_ ) {  ## if ERRFIL created but empty, return
	  unlink $ERRFIL;
      }
      else {  ##  there is an error message
	  unless ($suppress) {
	    $errormsg = $_;
            unlink $ERRFIL;
            }
      }  ## else there's an error message
    } ##  else there's an ERRFIL created

    if ($? || ($errormsg)) {  ## call err_routine if either $? or $errormsg
	if ($err_routine) {
	    &$err_routine($err_string,$err_log);
	}
	unshift @result, 1, "Problem running:\n $command:\n\n$errormsg\n" if ($errormsg);
	unshift @result, $?, "Problem running:\n $command\n" if (!$errormsg);
    } ## if any errors
    else {
	unshift @result, $?, "$command\n";
    }
    return @result;

} ## sub call

sub warn {

    local($string) = $_[0];
    local($log) = $_[1];
    $CURHANDLE = select(STDERR);
    print "-------------------------------------------------------------\n";
    print "WARNING:  $string\n";
    print "See full log in $log\n" if $log;
    print "-------------------------------------------------------------\n";
    select($CURHANDLE);
}

sub fdmp {
  local($infile) = $_[0];
  local($columns) = $_[1];
  return `fdump prhead=no infile=${infile} outfile=STDOUT columns=\"${columns}\" showcol=no showrow=no showunit=no rows=- more=yes wrap=yes pagewidth=200`;
}

# far fewer contortions to use the wonderful sumpha :-)
sub mergepha {
  local($inlist) = $_[0];
  local($outfile) = $_[1];
  local($log) = $_[2];
  local($com,@mergelist,@out);

  @mergelist = split(' ',$inlist);
  # better use '@filename' syntax in case of 
  # really long names and/or lots of files
  open(FLST,">aux/sumphalist") || die "\nCannot open aux/filelist to write:  $!\nRex quitting.\n";
  foreach $file (@mergelist){print FLST "$file\n"}
  close(FLST);
  $com = "sumpha filelist=\@aux/sumphalist outfile=$outfile clobber=yes";
  @out = &call($com,"warn","Problem running: \n$com\n","$log");
  unlink("aux/sumphalist");
  return;
}


# contortions to use the annoying mathpha on up to 100 files
sub mergepha_old {
  local($inlist) = $_[0];
  local($outfile) = $_[1];
  local($log) = $_[2];
  local($com,$expr,$bigexpr,$n,$m,@out,@mergelist);
  
  @mergelist = split(' ',$inlist);
  $m = 0; # count middle level of merged files
  while ($m < (scalar(@mergelist)/10) ) {
    $n = 0;  # count bottom level of original spectra 
    $expr = "";
    # add ten spectra at a time into ${m}.pha 
    while (($n < 10) && ($n < (scalar(@mergelist) - $m*10 ) )) {
      $expr .= "$mergelist[$n+10*$m] + "; 
      $n++;
    }
    $expr =~ s/ \+ $//g;  # remove trailing " + "
    $com = "mathpha expr=\"$expr\" units=C outfil=r${m}.pha exposure=CALC areascal=\% ncomments=0 clobber=yes";
    
    @out = &call($com,"warn","Problem running: \n$com\n","$log");
    if (shift(@out)) { # if the merge failed
      open(MLOG,">>$log") or die "Cannot open $log to write: $!\n";
      print MLOG @out;
      close(MLOG);
    }
    else { # start making sencond round expression
      $bigexpr .= "r${m}.pha + ";
    }
    $m++;
  }
  # if there were more than 10 originally, add up all temp spectra 
  # (assumed to be fewer than 10;  no obsid has > 100 data files
  if ($m > 1) {
    $bigexpr =~ s/ \+ $//g;  # remove trailing " + "

    $com = "mathpha expr=\"$bigexpr\" units=C outfil=$outfile exposure=CALC areascal=\% ncomments=0 clobber=yes";
    @out = &call($com,"warn","Problem running: \n$com\n","$log");
    if (shift(@out)) {
      open(MLOG,">>$log") or die "Cannot open $log to write: $!\n";
      print MLOG @out;
      close(MLOG);
    }
    foreach ($i=0; $i < (scalar(@mergelist)/10); $i++) {
      unlink "r${i}.pha";
    }
  }
  else { # if there were fewer than 10 originally, you're done
    system("mv r0.pha $outfile");
  }
  system("cphead infile=$mergelist[0]+1 outfile=${outfile}+1 keyfil=aux/keyfile.txt"); 
  return;
}

sub mergelcs {
  local($inlist) = $_[0];
  local($outfile) = $_[1];
  local($log) = $_[2];
  local($com,@out);
  @mergelist = split(' ',$inlist);
  unlink("lc.list");
  open(LIST,">lc.list") or die "Cannot open lc.list to write:  $!\n";
  foreach (@mergelist) {
    print LIST "$_ \n";
  }
  close(LIST);
  $com = "fmerge lastkey=TSTOPI infiles=\@lc.list outfile=${outfile} columns=- clobber=yes ";
  @out = &call("$com","warn","Problem running: \n $com\n","$log");
  unlink("lc.list");
  return @out;
}