#! /bin/sh
# This is the LHEA perl script: /cvmfs/extras-fp7.egi.eu/extras/heasoft/swift/x86_64-unknown-linux-gnu-libc2.19-0/bin/xrtgrblc
# 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/swift/x86_64-unknown-linux-gnu-libc2.19-0/bin/xrtgrblc."
  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/swift/x86_64-unknown-linux-gnu-libc2.19-0/bin/xrtgrblc."
  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!
#-------------------------------------------------------------------------------
#!/usr/bin/perl -w
#
# xrtgrblc - 
#       
#       extracts an XRT light curve and fitted spectra from a series 
#       of Swift GRB observations. It does so using a wide variety of 
#       HEASoft/Swift tasks, and relatively few inputs. It was designed
#       to automate the somewhat cumbersome task of extracting a corrected, 
#       flux converted light curve from a potentially large set of Swift 
#       observations.
#
#
# Modification History:
#
# 04/25/2008 - 1.0 - CAP
#              initial version
#
# 10/03/2008 - 1.1 - CAP - 
#              Many changes:
#
#               Bug Fixes:
#                * Upper limit handling and removal of low count bins
#                * Minor platform dependent bug in writing FITS info table
#                * Trap 0.0 for normalization in Xspec fits
#                * Trap various possible divides by zero
#                * Fixed bug in background PHA calculation
#                * Fixed empty list bug in PHA creation
#                * Trap failure to produce exposure map for a given interval
#                * Fixed double correction in flux light curves
#                * Fixed high memory usage when 'wtframetime' is small
#
#               Modifications:
#                * Changed Xspec PSF/Pileup correction to include expomap in ARF
#                * Added new method for correction:
#                   now can use XSPEC || XRTLCCORR
#                * Changed the behaviour when 'clobbering'
#                * New parameter 'cleanspec'
#                   == 'no': intermediate events, spectra, ARFS, and expo maps
#                            are not cleaned 
#                * Merged output into single FITS for each Mode/Energy range
#                  combination
#                * Turned off centroiding by default
#
# 12/01/2008 - 1.2 - CAP -
#              Bug Fixes:
#               * Binning code was not properly merging bins
#               * Trap zero-rows errors in fstatistic, and Xspec
#               * Trap zero counts in Ximage
#              
#              Enhancements/Modifications:
#               * Changed default PSF correction method to XRTLCCORR
#               * Increased background region size for PC mode
#               * Increased pile-up region for very high WT count rates.
#               * Added code to shift PC background regions, if they would go off-chip
#               * New hidden parameter 'object', defines OBJECT keyword in output files
#                 object='DEFAULT' uses first OBJECT keyword found
#               * New hidden parameter 'usetrigtime', default = 'yes'
#                  if ='yes', light curves are offset by TRIGTIME keyword or
#                  trigtime parameter
#                  if ='no', light curves are offset by earliest TSTART in input files
#               * New hidden parameter 'usetargid', default = 'yes'
#                 determines if TARGID keyword is written to output files
#
# 03/01/2009 - 1.3 - CAP -
#
#              Bug Fixes:
#               * Use maketime expression (param maketimeexpr) to generate a full
#                 GTI for all PC data. This total GTI is used to extract the
#                 total PC events used for detecting field sources.
#               * Changed PHA handling to not ignore PHAs with zero counts. This
#                 effects the normalized count rate in Xspec, since the exposure
#                 time was not correct.
#               * Fixed bug in updating orbit stop times after determining where to
#                 split orbits.
#               * Fixed wrong calculation of split energy band flux conversion -
#                 now uses model flux in split energy bands divided by
#                 corresponding count rate
#               * Work-around for file limits in mathpha
#               * Work-around mgtime file limits.
#               * Fixed a handful of undeclared variables
#              
#              Modifications:
#               * Lowered distance at which a detected source is considered a field
#                 source from .01 degrees to 0.0028 degrees.
#               * Update DATE-OBS/DATE-END keywords
#
#
# 06/29/2009 - 1.4 - CAP -
#
#              Bug Fixes:
#               * Fixed field source distance calculation, and region size calculation
#               * Fixed bug when calculating ARFs for intervals with zero counts.
#              
#              Modifications:
#               * Added fldsrcdist parameter to let the user decide how close
#                 to the input coordinates to allow field sources to be considered
#                 field sources. Inside of this parameter, any source detected
#                 is considered to be the GRB.
#               * Added cutlowbins(wt|pc) parameters to allow finer control
#                 over bin cutting
#               * Added trigfrombat parameter - used to indicate whether the 
#                 input trigger time is from the BAT or some other observatory.
#                 Only really affects plots, and header keyword comment for TRIGTIME
#
# 1/12/2012 - 1.6 - T. A. Reichard
#              Bug fixes:
#               * Forces Xspec to exit instead of wait for keyboard input when 
#                 a bad fit (too few good bins compared to the number of fit 
#                 parameters) occurs.
#               * Fixed a type-casting error in the PA_PNT (roll angle) 
#                 keyword value to avoid abrupt termination.
# 4/10/2014 - 1.8 - T. A. Reichard
#              New features:
#               * Added the lccorravg parameter, which can toggle between 
#                 using a mean PSF correction (pre-existing default) or
#                 a time-dependent PSF correction (new). The latter option
#                 runs more slowly but may be more accurate for periods when
#                 the attitude (and thus also the PSF correction) has not 
#                 stabilized.
#               * Added the pcreglist and wtreglist parameters, which can
#                 toggle between old (List #1) and new (List #2) source
#                 and background extraction regions in each mode.
#               * When light curve bins have a negative net rate, skip the 
#                 bin instead of keeping it. This change avoids negative
#                 correction factors.
#               * Added user-specified flux conversion factor parameters 
#                 6 (pc|wt)(hi|lo|)fluxconv to allow a user to specify
#                 the factors instead of using those generated by Xspec fits.

use strict;

#
# global variables
#
my @clean = ( );      # global list of files to clean if asked to
my @cleanspec = ( );  # directories for each obsid to clean

my $taskName    = 'xrtgrblc';
my $taskVersion = '1.8';

# instant gratification
$| = 1;

# input parameters
my %params = ( 
    'evtfiles'      => [],    # clean event file list
    'mkffiles'      => [],    # makefilter files
    'xhdfiles'      => [],    # XRT header housekeeping files
    'attfiles'      => [],    # attitude files
    'outdir'        => '',    # output directory
    'outstem'       => '',    # rootname of all output files
    'object'        => '',    # object name
    'srcra'         => '',    # source RA
    'srcdec'        => '',    # source Dec
    'fldsrcdist'    => 0.0,   # min dist to consider a source a field one
    'usecentroid'   => 0,     # refine position with xrtcentroid?
    'centroidhbox'  => 0.4,   # half-box length for xrtcentroid
    'psfcorrmeth'   => '',    # method for psf correction
    'lccorravg'     => 1,     # use average (1) or time-dependent (0) xrtlccorr PSF correction
    'maketimeexpr'  => '',    # GTI expression for maketime 
    'detthresh'     => 0.0,   # sigma > noise considered a detection
    'minpccnts'     => 0.0,   # minimum PC mode source counts to allow
    'minwtcnts'     => 0.0,   # minimum WT mode source counts to allow
    'pcreglist'    => 1,     # PC source/background regions table number
    'wtreglist'    => 1,     # WT source/background regions table number
    'wtframetime'   => 0.0,   # custom WT mode "unbinned" bin duration
    'pcframetime'   => 0.0,   # custom PC mode "unbinned" bin duration
    'minenergy'     => 0.0,   # minimum energy to consider
    'maxenergy'     => 0.0,   # maximum energy to consider
    'splitenergy'   => 0.0,   # where to split energy bands
    'splitorbits'   => 0,     # split orbits into sub-orbit intervals
    'minredchisq'   => 0.0,   # minimum red. chi-squared
    'maxredchisq'   => 0.0,   # maximum red. chi-squared
    'pcfluxconv'    => 0.0,   # PC flux conversion factor for total energy band
    'pclofluxconv'  => 0.0,   # PC flux conversion factor for lowest energy band
    'pchifluxconv'  => 0.0,   # PC flux conversion factor for highest energy band
    'wtfluxconv'    => 0.0,   # WT flux conversion factor for total energy band
    'wtlofluxconv'  => 0.0,   # WT flux conversion factor for lowest energy band
    'wthifluxconv'  => 0.0,   # WT flux conversion factor for highest energy band
    'usetrigtime'   => 0,     # use trig time as time reference?
    'trigtime'      => 0.0,   # GRB trigger time
    'trigfrombat'   => 0,     # GRB trigger time is from the BAT?
    'usetargid'     => 0,     # use target ID in output files and plots?
    'nhmap'         => 0,     # which map to use in nh (0=LAB, 1=DL, -1 user)
    'usernh'        => 0.0,   # user specified NH value (cm^-2)
    'chatter'       => 0,     # chatinness level
    'bincurve'      => 0,     # create binned curve also?
    'bintype'       => 0,     # type of binning
    'minbincts'     => 0.0,   # minimum counts/bin in binned l.c.
    'pcsteprates'   => '',    # list of "upper rate boundary:cts/bin"
    'wtsteprates'   => '',    # list of "upper rate boundary:cts/bin"
    'cutlowbins'    => 0,     # drop bins with < minbincts counts in them?
    'cutlowbinspc'  => 0,     # cut low bins (PC mode only)
    'cutlowbinswt'  => 0,     # cut low bins (WT mode only)
    'gapintv'       => 0.0,   # minimum time in seconds considered a gap
    'plotcurves'    => 0,     # generate qdp files and plots
    'plotftype'     => '',    # plot type (valid PGPLOT file device)
    'plotfext'      => '',    # NOT A PARAM, but is the plot file extension
    'clean'         => 0,     # clean un-needed output files
    'cleanspec'     => 0,     # clean-up evt/spectra for each interval?
    'clobber'       => 0,     # clobber output dir?
    'history'       => 0      # add history?
);

# libs
use HEACORE::HEAUTILS;

use Astro::FITS::CFITSIO qw( :longnames :constants );

use File::Spec::Functions qw( :ALL );
use File::Basename;
use File::Path;
use File::Copy;
use Cwd;
use Math::Trig;

use HEACORE::PIL;

use libswxrtgrblc;
    
#
# HEADAS wrapper
#
use HEACORE::HEAINIT;
exit headas_main( \&xrtgrblc );

#
# the xrtgrblc subroutine
#
sub xrtgrblc {
    
    # task name and version
    set_toolname( $taskName );
    set_toolversion( $taskVersion );
    libswxrtgrblc::setprefix( "${taskName}_${taskVersion}: " );
    libswxrtgrblc::setchat( 2 );
    
    # status variable
    my $stat = 0;
    
    # get and check all input parameters
    $stat = getParams( \%params );
    quit( $stat ) unless $stat == 0;
  
    # parse the steprates parameters
    my %steprates;
    ( $stat, %steprates ) = parseStepRates( );
    quit( $stat ) unless $stat == 0;

    # nice little startup message
    my $time  = localtime( time );
    my $intro = "---- Starting $taskName v$taskVersion at $time ----";
    my $dashs = "-" x length( $intro );
    chat( 1, "$dashs\n$intro\n$dashs\n\n" );
    
    # clobber
    if ( -d $params{outdir} && !$params{clobber} ) {
        error( 1, "output dir $params{outdir} exists! Check clobber param\n" );
        quit( -1 );
    } else {
        eval { mkpath( $params{outdir} ); };
        if ( $@ ) {
            error( 1, "failed to make output dir $params{outdir}\n" );
            quit( -1 );
        }
    }

    # check the plot device if plotting
    if ( $params{plotcurves} ) {
        $stat = getPlotDevice( );
        quit( $stat ) unless $stat == 0;
    }
    
    # PI for XRT is roughly 100*Energy, and "approximately linear"
    my $piMin = 100.0 * $params{minenergy};
    my $piMax = 100.0 * $params{maxenergy};
    my $piBrk = 100.0 * $params{splitenergy};

    # get the source nH and scale for xspec
    my $nh;
    if ( $params{nhmap} >= 0 ) {
        ( $stat, $nh ) = getNH( $params{srcra}, $params{srcdec}, 2000, $params{nhmap} );
        quit( $stat ) unless $stat == 0;
    } else {
        $nh = $params{usernh};
    }
    
    my $xspecNH = $nh / 1.0e22;

    # check the FOV of each event file
    $stat = checkEventsFOV( \@{$params{evtfiles}}, $params{srcra}, $params{srcdec} );
    quit( $stat ) unless $stat == 0;

    # everything, essentially, goes in here
    my %evtData = ( );

    # list of pc event files only
    my @pcEvtsList = ( );

    # list of GTIs for combining PC events
    my @pcGTIs = ( );
    
    # setup the output dirs
    for ( my $i = 0; $i < @{$params{evtfiles}}; $i++ ) {

        my $obsDir   = undef;    # where stuff for this event file goes
        my $obsid    = undef;    # OBS_ID for this event file
        my $segnum   = undef;    # segment number
        my $mode     = undef;    # PC or WT
        my $submode  = undef;    # W1, W2, W3, W4 or W5
        my $midpix   = undef;    # midpoint of the PC mode window
        my $pa_pnt   = undef;    # Swift roll angle
        my $startd   = undef;    # DATE-OBS
        my $stopd    = undef;    # DATE-END
        my $framtime = undef;    # Frame time
        my $bint     = undef;    # "unbinned" l.c. time
        my $lcthresh = undef;   # minimum FRACEXP to not lose any events
        my $grade    = undef;    # grade selection string
        my $vsub     = undef;    # XRTVSUB

        # make the output directory (getting obsid/mode/submode in the process)
        $stat = createOutputDir( $params{evtfiles}[ $i ], $params{outdir}, \$obsDir, 
                                 \$startd, \$stopd, \$obsid, \$mode, \$submode, 
                                 \$midpix, \$pa_pnt, \$framtime, \$grade, \$vsub );
        quit( $stat ) unless $stat == 0;
        push @cleanspec, $obsDir;

        # fix the bin time based on input param, but don't let bin duration
        # go below frame time
        if ( uc( $mode ) eq 'WT' && $params{wtframetime} >= $framtime ) {
            $bint = $params{wtframetime};
        } elsif ( uc( $mode ) eq 'PC' && $params{pcframetime} >= $framtime ) {
            $bint = $params{pcframetime};
        } else {
            $bint = $framtime;
        }

        # hash each mode on full path of filename
        if ( exists $evtData{$params{evtfiles}[ $i ]} ) {
            error( 1, "duplicate event file $params{evtfiles}[ $i ] for $obsid\n" );
            quit( -1 );
        }

        # parse the filename
        my $rn;
        Astro::FITS::CFITSIO::fits_parse_rootname( $params{evtfiles}[ $i ], $rn, $stat );
        quit( $stat ) unless $stat == 0;

        my $bn = basename( $rn );
        $bn    =~ s/^(.*?)(\.gz)?$/$1/;

        # get the segment number
        $segnum = $obsid;
        $segnum =~ s/['"]//g;
        $segnum =~ s/^\d{8}//;
        
        # get relevant CALDB files
        my $teldef = "";
        my $align  = "";
        my $rmf    = "";
        
        $stat = getCalDBFiles( $mode, $startd, $stopd, $grade, $vsub, 
                               \$teldef, \$align, \$rmf ); 
        quit( $stat ) unless $stat == 0;
    
        
        # per event file bookeeping data (written to fits info file later)
        $evtData{$params{evtfiles}[ $i ]} = { 
            'evtsfile'   => $params{evtfiles}[ $i ],
            'base'       => $bn,
            'dir'        => $obsDir,
            'stem'       => "${segnum}_${mode}${submode}",
            'mkffile'    => $params{mkffiles}[ $i ],
            'xhdfile'    => $params{xhdfiles}[ $i ],
            'attfile'    => $params{attfiles}[ $i ],
            'mode'       => $mode,
            'submode'    => $submode,
            'bintime'    => $bint,
            'framtime'   => $framtime,
            'midpix'     => $midpix,
            'pa_pnt'     => $pa_pnt,
            'obsid'      => $obsid,
            'teldef'     => $teldef,
            'align'      => $align,
            'rmf'        => $rmf,
            'lcthresh'   => $lcthresh,
            'start'      => [ ],                # start time of each interval
            'stop'       => [ ],                # stop time of each interval
            'orbit'      => [ ],                # actual orbit number
            'pa_pnts'    => [ ],                # roll angles (WT only) per-orbit
            'bkgreg'     => [ ],                # background regions
            'srcreg'     => [ ],                # source regions
            'orbevts'    => [ ],                # per-orbit events file (for expomap)
            'expomap'    => [ ],                # per-orbit exposure map
            'srccnts'    => [ ],                # source counts
            'hsrccnts'   => [ ],                # source counts hard band
            'ssrccnts'   => [ ],                # source counts soft band
            'uplimit'    => [ ],                # source uplimit
            'huplimit'   => [ ],                # source uplimit hard band
            'suplimit'   => [ ],                # source uplimit soft band
            'bkgcnts'    => [ ],                # bkg counts
            'hbkgcnts'   => [ ],                # bkg counts hard band
            'sbkgcnts'   => [ ],                # bkg counts soft band
            'engconv'    => [ ],                # energy conversion
            'sengconv'   => [ ],                # soft energy conversion
            'hengconv'   => [ ],                # hard energy conversion
            'psfcorr'    => [ ],                # mean psf area correction
            'psfcorrfacts' => [ ],              # array of psf area correction factors
            'psfcorrtimes' => [ ],              # array of psf area correction times
            'psfcorr'    => [ ],                # psf area correction
            'expcorr'    => [ ],                # exposure correction
            'backscl'    => [ ],                # background area scale
            'backcorr'   => [ ],                # total bkg correction
            'exposure'   => [ ],                # orbit exposure time
            'detected'   => [ ],                # source detected?
            'chisq'      => [ ],                # chi-squared of P.L. fit
            'dof'        => [ ],                # deg. of free. of P.L. fit
            'plindex'    => [ ]                 # power-law index
        };

        # run make time to get orbits (as defined by maketimeexpr param)
        my $gtiFile = undef;

        $stat = runMakeTime( $evtData{$params{evtfiles}[ $i ]}, $params{maketimeexpr}, \$gtiFile );
        quit( $stat ) unless $stat == 0;

        # save for later
        $evtData{$params{evtfiles}[ $i ]}->{gtifile} = $gtiFile;

        # if PC mode, add to list of GTIs for combined event file
        if ( $mode eq 'pc' ) {
            push @pcGTIs, $gtiFile;
	    push @pcEvtsList, $params{evtfiles}[ $i ];
        }
    }
   
    # check various timing keywords in each event file
    ( $stat, $params{trigtime} ) = checkTiming( \%evtData, $params{trigtime} );
    quit( $stat ) unless $stat == 0;

    # generate a single events file with all PC mode and detect all sources
    # refine ra/dec if applicable
    my $pcEvts;
    my @detSources = ( );       # detected sources for PC mode
    my @fldRegions = ( );       # field source regions for PC mode

    if ( @pcEvtsList ) {

        # merge the PC mode GTIs
        my $pcGTI = catfile( $params{outdir}, "$params{outstem}_pcgti.fits" );
        push @clean, $pcGTI;

        $stat = orGTIs( $pcGTI, \@pcGTIs );
        quit( $stat ) unless $stat == 0;
   
        # combine events, including only our orbit selection
        ( $stat, $pcEvts ) = combinePCEventsFiles( \@pcEvtsList, $pcGTI, $params{outdir}, 
                                                   $params{outstem}, $piMin, $piMax );
        
        quit( $stat ) unless $stat == 0;
       
        # detect
        my $pcDetBase = basename( $pcEvts );
        $pcDetBase =~ s/\.evt$//;
        my $pcDetFile;
        ( $stat, $pcDetFile ) = detectSources( $pcEvts, $params{outdir}, $pcDetBase,
                                               $params{srcra}, $params{srcdec}, 
                                               $piMin, $piMax );
        
        quit( $stat ) unless $stat == 0;                     
        
        # centroid if asked
        if ( $params{usecentroid} ) {
            my ( $nra, $ndec );
            ( $stat, $nra, $ndec ) = runXRTCentroid( $pcEvts, $pcDetBase, $params{outdir},
                                                     $params{srcra}, $params{srcdec}, 
                                                     $params{centroidhbox} );
            quit( $stat ) unless $stat == 0;
                
            $params{srcra}  = $nra;
            $params{srcdec} = $ndec;
        }
       
        # read detected sources if file exists, otherwise assume GRB is a source
        if ( -e $pcDetFile ) {
            $stat = readXImageDetectFile( $pcDetFile, \@detSources );
            quit( $stat ) unless $stat == 0;
        }

        # in case we don't detect the GRB, put it in the source list
        my $det = {
            RA      => $params{srcra},
            DEC     => $params{srcdec}
        };
        push @detSources, $det;

    }

    # get extraction regions
    foreach my $key ( keys %evtData ) {

        my $evt = $evtData{$key};

        ## run make time to get orbits (as defined by maketimeexpr param)
        #my $gtiFile = undef;

        #$stat = runMakeTime( $evt, $params{maketimeexpr}, \$gtiFile );
        #quit( $stat ) unless $stat == 0;
        
        # get the orbit start and stop times
        $stat = getOrbits( $evt->{gtifile}, \@{$evt->{start}}, \@{$evt->{stop}} );
        quit( $stat ) unless $stat == 0;
        
        # determine the per orbit extraction regions
        if ( $evt->{mode} eq 'pc' ) {
            $stat = getExtractionRegionsPC( $evt, $params{pcreglist}, 
					    $params{srcra}, 
                                            $params{srcdec}, $piMin, $piMax,
                                            \@detSources, $params{splitorbits} );
            quit( $stat ) unless $stat == 0;
            
        # WT
        } else {
            $stat = getExtractionRegionsWT( $evt, $params{wtreglist}, 
					    $params{srcra}, 
                                            $params{srcdec}, $piMin, 
                                            $piMax, $params{splitorbits} );
            quit( $stat ) unless $stat == 0;
        }
       
        # purge event files with no usable data
        if ( !@{$evt->{start}} ) {
            warnhi( 1, "no good intervals found in $evt->{evtsfile}\n" );
            delete $evtData{$key};
        }
    }
    
    # pop the GRB RA/Dec from the detected sources list
    pop @detSources;

    # get a list of exclusion regions - only those detected sources
    # that are within the maximum background annulus (310" currently)
    foreach my $detsource ( @detSources ) {
        if ( exists $detsource->{FLDDIST} && $detsource->{FLDDIST} <= 310.0 ) {
            push @fldRegions, $detsource->{REGION};
        }
    }

    # check if we have anything good to do
    if ( !keys %evtData ) {
        error( 1, "no good data found\n" );
        quit( -1 );
    }
   
    ##########################################
    # below is the real light curve extraction
    ##########################################

    # info file
    my $infoFile = catfile( $params{outdir}, "$params{outstem}_info.fits" );

    # list of GTIs to merge
    my @gtis = ( );
    
    # arf/pha lists for merging
    my %arfs    = ( );
    my %exps    = ( );
    my %phas    = ( );
    my %bkgphas = ( );
    my %bkgscls = ( );
    my %bkgexps = ( );
    my %totexpo = ( );
    my %totbkgexpo = ( );
   
    # start with full energy range
    my $emin = $params{minenergy};
    my $emax = $params{maxenergy};
    my $pmin = $piMin;
    my $pmax = $piMax;

    # ratio plot min/max
    my $ratiomin =  1.0e20;
    my $ratiomax = -1.0e20;

    # plot files
    my $lcplotfile  = catfile( $params{outdir}, "$params{outstem}_xpwetsrab.qdp" );
    headas_clobberfile( $lcplotfile );
    if ( -e $lcplotfile ) {
        error( 1, "failed to clobber file $lcplotfile\n" );
        quit( -1 );
    }
    my $ratplotfile = catfile( $params{outdir}, "$params{outstem}_xpwetsrrb.qdp" );
    headas_clobberfile( $ratplotfile );
    if ( -e $ratplotfile ) {
        error( 1, "failed to clobber file $ratplotfile\n" );
        quit( -1 );
    }

    # files that need keyword fixing
    my %fixFits = ( );
    
    # get spectrum corrections only from full band
    my $band = 'fb';

    # do the modes separately
    my $mode = 'wt';

    # from here may be repeated depending on the band splitting and modes used
EXTRACT:

    # conversions/corrections from previous orbit/event list
    my $lastEconv  = 5.0e-11;   # default energy conversion, probably close
    my $lastSEconv = 5.0e-11;   # default soft-band energy conversion
    my $lastHEconv = 5.0e-11;   # default hard-band energy conversion

    my $userflux = 0;

    if ( $params{"${mode}fluxconv"} > 0.0 ) {
        $userflux = 1;
    }
    if ( $userflux && $params{splitenergy} > 0.0 &&
         $params{"${mode}hifluxconv"} <= 0 ) {
        $userflux = 0;
    }
    if ( $userflux && $params{splitenergy} > 0.0 &&
         $params{"${mode}lofluxconv"} <= 0 ) {
        $userflux = 0;
    }
    if ( $userflux ) {
        $lastEconv = $params{"${mode}fluxconv"};
        $lastSEconv = $params{"${mode}lofluxconv"};
        $lastHEconv = $params{"${mode}hifluxconv"};
    }


    # set mode- and band-specific filename base
    my $modebase = "$params{outstem}_x${mode}";
    my $bandbase;
    if ( $band eq 'fb' ) {
        $bandbase = sprintf "${modebase}et";
    } elsif ( $band eq 'soft' ) {
        $bandbase = sprintf "${modebase}e1";
    } elsif ( $band eq 'hard' ) {
        $bandbase = sprintf "${modebase}e2";
    }
    
    # various light curve files for each band
    my $lcbase    = catfile( $params{outdir}, $bandbase );

    my $allLCFile = $lcbase . "sra.lc";

    # binned versions
    my $allBLCFile = $lcbase . "srab.lc";
    push @clean, $allBLCFile;
   
    # the gti
    my $allGTIFile = $lcbase . ".gti";
    push @clean, $allGTIFile;

    # clear the GTI list
    @gtis = ( );

    # hardness ratio (band is listed as et)
    my $ratBLCFile = catfile( $params{outdir}, "${modebase}etsrrb.lc" );
    my $hardLCFile;

    # clear arf/pha lists for merging
    %arfs    = ( );
    %exps    = ( );
    %phas    = ( );
    %bkgphas = ( );
    %bkgscls = ( );
    %bkgexps = ( );
    %totexpo = ( );
    %totbkgexpo = ( );
   
    # corresponding filenames - with muck for mathpha
    foreach my $clf ( $allLCFile, $allBLCFile, $allGTIFile, $ratBLCFile ) {
        headas_clobberfile( $clf );
        if ( -e $clf ) {
            error( 1, "failed to clobber file $clf\n" );
            quit( -1 );
        }
    }

    # template light curve (for various instrumental keywords)
    my $lctempl;
    my $rattempl;

    # get empty light curve
    my %rawLC = getEmptyLC( );

    # binned curves
    my ( $rawBLC, $hardBLC );

    # FITS pointer that gets passed around a little
    my $allfptr;

    # loop over each data set, in time sorted order
    # at least, time sorted with respect to start time
    # of first GTI interval - may not always work,
    # so the final curve is time sorted.
    my $first = 1;
    foreach my $key ( sort {
                             $evtData{$a}{start}[0] <=> 
                             $evtData{$b}{start}[0]
                           } 
                      keys %evtData ) {

        # we do some splicing in certain cases, so check
        # for existence of this hash value
        if ( !exists $evtData{$key} ) { next; }
        my $evt = $evtData{$key};
       
        # no mode mixing
        if ( $evt->{mode} ne $mode ) { next; }

        # do each interval
        for ( my $i = 0; $i < @{$evt->{start}}; $i++ ) {
           
            my ( $cmd, $out );
            
            # file labels
            my $orb = sprintf( "%02d", $i + 1 );

            my $grp;
            if ( $band eq 'fb' ) {
                $grp = $orb;
            } else {
                $grp = sprintf( "%02d_$band", $i + 1 );
            }

            # cat the source/bkg regions with negative field source regions
            my $srcreg;
            my $bkgreg;
            if ( $evt->{mode} eq 'pc' ) {
                $srcreg = join "\n", $evt->{srcreg}[ $i ], @fldRegions;
                $bkgreg = join "\n", $evt->{bkgreg}[ $i ], @fldRegions;
            } else {
                $srcreg = $evt->{srcreg}[ $i ];
                $bkgreg = $evt->{bkgreg}[ $i ];
            }

            # temp filenames - both per interval and per band/interval
            my ( $srclc, $bkglc, $srcspec, $bkgspec, $orbevts, $orbbkgevts, $orbexpo, $orbinst );
            my $orbbase = catfile( $evt->{dir}, "$evt->{stem}" );
            
            # source light curve
            $srclc   = $orbbase . "_${grp}.lc";
            
            # background light curve
            $bkglc   = $orbbase . "_${grp}_bkg.lc";
            
            # exposure and instrument map for this interval
            $orbexpo = $orbbase . "_${orb}_ex.img";
            $orbinst = $orbbase . "_${orb}_rawinstr.img.gz";
            
            push @cleanspec, $srclc, $bkglc, $orbexpo, $orbinst;
            
            # spectra and event files per interval
            $orbevts    = $orbbase . "_${orb}.evt";
            $orbbkgevts = $orbbase . "_${orb}_bkg.evt";
            if ( $band eq 'fb' ) {
                $srcspec    = $orbbase . "_${orb}.pha";
                $bkgspec    = $orbbase . "_${orb}_bkg.pha";
                headas_clobberfile( $orbevts );
                headas_clobberfile( $orbbkgevts );
                if ( -e $orbevts || -e $orbbkgevts ) {
                    error( 1, "failed to clobber files\n" );
                    quit( -1 );
                }
                push @cleanspec, $srcspec, $bkgspec, $orbevts, $orbbkgevts;
            } else {
                $srcspec = $bkgspec = undef;
            }
            foreach my $clf ( $srclc, $bkglc, $srcspec, $bkgspec ) {
                defined $clf && headas_clobberfile( $clf );
                if ( defined $clf && -e $clf ) {
                    error( 1, "failed to clobber file $clf\n" );
                    quit( -1 );
                }
            }
            
            # extract the real light curve from the defined region - source and bkg
            # extract a spectrum, events also, if this is the full band
            # and save the regions to disk
            chat( 2, "extracting source and bkg light curve and spectrum from ",
                  "$evt->{base}\n" );
            chat( 2, "from SCC times $evt->{start}[ $i ] to ",
                  "$evt->{stop}[ $i ]\n" );
            chat( 2, "with energy range ", sprintf( "%.2f", $emin ), "keV to ",
                  sprintf( "%.2f", $emax ), "keV\n" );
             
            # reuse previously extracted event file if extant
            my ( $srccnts, $srcarea, $bkgcnts, $bkgarea );
            my $useevents;
            if ( $band eq 'fb' ) {
                $useevents = $evt->{evtsfile};
            } else {
                $useevents = $orbevts;
                $orbevts = undef;
            }

            ( $stat, $srccnts ) = extractLC( $useevents, "$evt->{stem}_${orb}", 
                                             $evt->{dir}, $pmin, $pmax, 
                                             $evt->{start}[ $i ], $evt->{stop}[ $i ], 
                                             0.0, $srcreg, 1, $evt->{bintime},
                                             $srclc, $srcspec, $orbevts );
            quit( $stat ) unless $stat == 0;

            if ( $band eq 'fb' ) {
                $useevents = $evt->{evtsfile};
            } else {
                $useevents = $orbbkgevts;
                $orbbkgevts = undef;
            }
            ( $stat, $bkgcnts ) = extractLC( $useevents, "$evt->{stem}_${orb}_bkg",
                                             $evt->{dir}, $pmin, $pmax, 
                                             $evt->{start}[ $i ], $evt->{stop}[ $i ], 
                                             0.0, $bkgreg, 1, $evt->{bintime},
                                             $bkglc, $bkgspec, $orbbkgevts );
            quit( $stat ) unless $stat == 0;

            my $srcregFile = catfile( $evt->{dir}, "$evt->{stem}_${orb}.reg" );
            my $bkgregFile = catfile( $evt->{dir}, "$evt->{stem}_${orb}_bkg.reg" );
            push @cleanspec, $srcregFile, $bkgregFile;

            # get the area of source + bkg
            ( $stat, $srcarea ) = getFltFKey( "$srclc+1", 'NPIXSOU' );
            quit( $stat ) unless $stat == 0;
            
            ( $stat, $bkgarea ) = getFltFKey( "$bkglc+1", 'NPIXSOU' );
            quit( $stat ) unless $stat == 0;

            my $backscale;
            if ( $bkgarea > 0.0 ) {
                $backscale = $srcarea / $bkgarea;
            } else {
                $backscale = 0.0;
            }

            # For the full band, check detection and get conversions/corrections
            if ( $band eq 'fb' ) {
                my ( $detected, $dospecanal, $snr, $econv, $heconv, $seconv, 
                     $psfcorrfacts, $psfcorrtimes, $exposure, $expcorr, $bkgexpcorr, $plind, $dof, $chisq );

                # calc S/N ratio
                my $stot  = $srccnts - $bkgcnts * $backscale;
                my $denom = sqrt( $srccnts + $bkgcnts * $backscale**2.0 );
                if ( $denom == 0.0 ) {
                    $snr = 0.0;
                } else {
                    $snr = $stot / $denom;
                }

                # only consider a detection if S/N is high enough, and
                # only do spectral analysis if we have "enough counts"
                # (user determined)
                # the S/N is checked again during the binning phase
                $detected = ( $snr >= $params{detthresh} );
                if ( $srccnts < $params{"min${mode}cnts"} ) {
                    $dospecanal = 0;
                } else {
                    $dospecanal = $detected;
                }

                # generate exposure map if there are any source or bkg counts
                my %xrtexpomap = (
                    attfile       => $evt->{attfile},
                    hdfile        => $evt->{xhdfile},
                    stemout       => "$evt->{stem}_$orb",
                    fovfile       => 'CALDB',
                    vigfile       => 'CALDB',
                    vigflag       => 'no',
                    pcnframe      => 0,
                    wtnframe      => 0,
                    checkattitude => 'no',
                    outdir        => $evt->{dir},
                    teldef        => $evt->{teldef},
                    cleanup       => 'yes',
                    history       => 'yes',
                    clobber       => 'yes',
                    chatter       => $params{chatter}
                );
                if ( $srccnts > 0 ) {
                    $xrtexpomap{infile} = $orbevts;
                } else {
                    $xrtexpomap{infile} = $orbbkgevts;
                }

                $cmd = genCmd( 'xrtexpomap', \%xrtexpomap ); 
                chat( 2, "generating exposure map for $orbevts\n" );
                $stat = runSystem( $cmd );
                if ( $stat != 0 || !-e $orbexpo ) {
                    error( 1, "failed to produce exposure map for $orbevts\n",
                              "this orbit will be skipped\n" );
                    splice @{$evt->{start}}, $i, 1;
                    splice @{$evt->{stop}}, $i, 1;
                    splice @{$evt->{srcreg}}, $i, 1;
                    splice @{$evt->{bkgreg}}, $i, 1;
                    splice @{$evt->{orbit}}, $i, 1;
                    splice @{$evt->{pa_pnts}}, $i, 1;
                    if ( !@{$evt->{start}} ) {
                        delete $evtData{$key};
                    }
                    $i--;
                    $stat = 0;
                    next;
                }

                # make an arf with and without psf corrections if there are any counts
                my $specarfp = catfile( $evt->{dir}, "$evt->{stem}_${orb}_psf.arf" );
                my $specarfn = catfile( $evt->{dir}, "$evt->{stem}_${orb}_nopsf.arf" );
                push @cleanspec, $specarfp, $specarfn;

                if ( $srccnts > 0 ) {
                    $stat = makeArfs( $srcspec, $orbexpo, $specarfp, $specarfn );
                    quit( $stat ) unless $stat == 0;
                } else {
                    $dospecanal = 0;
                }
                    
                # get various conversions/corrections
                if ( $dospecanal ) {
		    
                    my $stats;
                    if ( uc( $params{psfcorrmeth} ) eq 'XSPEC' || !$userflux ) {
			# first group the spectrum to 20cts/bin (if possible)
			my $grpspec  = catfile( $evt->{dir}, "$evt->{stem}_${orb}_grp.pha" );

			push @cleanspec, $grpspec;
			
			$stat = groupSpectrum( $srcspec, $grpspec, $params{chatter} );
			
			# fit the spectra to get energy/psf conversions
			( $stat, $stats ) = 
			    getSpecCorrections( $evt->{dir}, "$evt->{stem}_$orb",
						$grpspec, $bkgspec, $specarfp, $specarfn, 
						$evt->{rmf}, $xspecNH, $params{minenergy},
						$params{maxenergy}, $params{splitenergy} );
			
			quit( $stat ) unless $stat == 0;
			$chisq   = $stats->{chisq};
			$dof     = $stats->{dof};
			$plind   = $stats->{plind};
                    } else {
                        $stats = {
                            econv   => $lastEconv,
                            heconv  => $lastHEconv,
                            seconv  => $lastSEconv,
                            psfcorr => 1.0,
                            plind   => 0.0,
                            plindlo => 0.0,
                            plindhi => 0.0,
                            dof     => 1,
                            chisq   => 1
                        };
                        $chisq = FLT_NULL;
                        $dof   = FLT_NULL;
                        $plind = FLT_NULL;
                    }

                    # use xrtlccorr to get PSF correction if asked to
                    if ( uc( $params{psfcorrmeth} ) eq 'XRTLCCORR' ) {
                        my $pa = $evt->{mode} eq 'pc' ? 1 : $evt->{pa_pnts}[ $i ];
                        ( $stat, $stats->{psfcorrfacts}, $stats->{psfcorrtimes} ) = 
                          getXrtlccorrCorrection( $mode, $srclc, $orbevts, $orbinst, $evt->{teldef},
                                                  $evt->{attfile}, $evt->{xhdfile}, $pa, $params{lccorravg} );
                          quit( $stat ) unless $stat == 0;
                    }

                    # check for fit failure
                    my $redchisq = $stats->{dof} == 0.0 ? FLT_NULL :
                                                    $stats->{chisq} / $stats->{dof};
                    if ( $stats->{econv} == FLT_NULL ||  
                         $redchisq > $params{maxredchisq} ||
                         $redchisq < $params{minredchisq} ) {
                        $econv   = $lastEconv;
                        $heconv  = $lastHEconv;
                        $seconv  = $lastSEconv;
                        
                        # fudge the PSF correction
                        if ( uc( $params{psfcorrmeth} ) ne 'XRTLCCORR' ) {
                            ( $stat, $psfcorrfacts, $psfcorrtimes ) = fudgePSFCorr( $specarfp, $specarfn, $orbevts );
                            quit( $stat ) unless $stat == 0;
                        } else {
                            $psfcorrfacts = $stats->{psfcorrfacts};
                            $psfcorrtimes = $stats->{psfcorrtimes};
                        }
                    } else {
                        $econv   = $stats->{econv};
                        $heconv  = $stats->{heconv};
                        $seconv  = $stats->{seconv};
                        $psfcorrfacts = $stats->{psfcorrfacts};
			$psfcorrtimes = $stats->{psfcorrtimes};
                    }

                # otherwise, if not detected, use the last conversion/correction
                } else {
                    $econv   = $lastEconv;
                    $heconv  = $lastHEconv;
                    $seconv  = $lastSEconv;
                    $chisq   = FLT_NULL;
                    $dof     = FLT_NULL;
                    $plind   = FLT_NULL;

                    # fudge the PSF correction
                    if ( $srccnts > 0 ) {

                        # use xrtlccorr to get PSF correction if asked to
                        if ( uc( $params{psfcorrmeth} ) eq 'XRTLCCORR' ) {
                            my $pa = $evt->{mode} eq 'pc' ? 1 : $evt->{pa_pnts}[ $i ];
                            ( $stat, $psfcorrfacts, $psfcorrtimes ) = 
                              getXrtlccorrCorrection( $mode, $srclc, $orbevts, $orbinst, $evt->{teldef},
                                                      $evt->{attfile}, $evt->{xhdfile}, $pa , $params{lccorravg});
                              quit( $stat ) unless $stat == 0;
                        } else {
                            ( $stat, $psfcorrfacts, $psfcorrtimes ) = fudgePSFCorr( $specarfp, $specarfn, $orbevts );
                            quit( $stat ) unless $stat == 0;
                        }
                    } else {
                        $psfcorrfacts = [1.0];
			$psfcorrtimes = [0.0];
                    }

                }

		my @sortedpsfcorrfacts = sort {$a <=> $b} @$psfcorrfacts;
                chat( 2,  "PSF correction for $srclc ranges from ",
		      sprintf( "%.3f", $sortedpsfcorrfacts[0] ),
		      " to ",
		      sprintf( "%.3f", $sortedpsfcorrfacts[-1] ) 
		    );
                chat( 2, "energy conversion for $srclc is: ", 
                      sprintf( "%.2e", $econv ), " erg/cm^2/count" );

                # get the nominal exposure time
                ( $stat, $exposure ) = getFltFKey( "$orbevts+1", 'EXPOSURE' );
                quit( $stat ) unless $stat == 0;

                # there is no exposure correction (it's already in the other corrections)
                $expcorr = $exposure;
                
                # add to the lists of PHAs
                # to merge if there are any source counts
                if ( -e $specarfp ) {

                    # save the arf and phas
                    # carefully, so we don't muck up mathpha with long filenames
                    my $fn = rel2abs( $srcspec );
                    $fn    = abs2rel( $fn, $params{outdir} );
                    if ( !exists $arfs{$evt->{rmf}} ) {
                        $arfs{$evt->{rmf}}  = [ $specarfp ];
                        $exps{$evt->{rmf}}  = [ $exposure ];
                        $phas{$evt->{rmf}}  = [ $fn ];
                        $totexpo{$evt->{rmf}} = $exposure;
                    } else {
                        push @{$arfs{$evt->{rmf}}}, $specarfp;
                        push @{$exps{$evt->{rmf}}}, $exposure;
                        push @{$phas{$evt->{rmf}}}, $fn;
                        $totexpo{$evt->{rmf}} += $exposure;
                    }

                    # save the pha and backscale if we have src cnts
                    $fn = rel2abs( $bkgspec );
                    $fn = abs2rel( $fn, $params{outdir} );
                    if ( !exists $bkgphas{$evt->{rmf}} ) {
                        $bkgphas{$evt->{rmf}}  = [ $fn ];
                        $bkgexps{$evt->{rmf}}  = [ $exposure ];
                        $bkgscls{$evt->{rmf}}  = [ $backscale ];
                        $totbkgexpo{$evt->{rmf}} = $exposure;
                    } else {
                        push @{$bkgphas{$evt->{rmf}}}, $fn;
                        push @{$bkgexps{$evt->{rmf}}}, $exposure;
                        push @{$bkgscls{$evt->{rmf}}}, $backscale;
                        $totbkgexpo{$evt->{rmf}} += $exposure;
                    }
                }

                # get the background exposure per pixel
                if ( -e $orbexpo ) {
                    ( $stat, $bkgexpcorr ) = getRegionCounts( $orbexpo, 
                                                              "$evt->{stem}_${orb}_bkg_expo",
                                                              $evt->{dir}, $bkgregFile );
                    quit( $stat ) unless $stat == 0;
                    if ( $bkgexpcorr <= 0.0 ) {
                        warnhi( 1, "no background exposure correction ",
                                   "required for interval $orb\n" );
                        $bkgexpcorr = $exposure;
                    }
                } else {
                    warnhi( 1, "no background exposure correction ",
                               "possible for interval $orb\n" );
                    $bkgexpcorr = $exposure;
                }

                # calc source exposure correction, and total background correction
                if ( $bkgexpcorr > 0.0 && $bkgarea > 0.0 ) {
                    $bkgexpcorr = ( $exposure * $srcarea ) / ( $bkgexpcorr * $bkgarea );
                } else {
                    $bkgexpcorr = 1.0;
                }
                if ( $expcorr > 0.0 ) {
                    $expcorr = $exposure / $expcorr;
                } else {
                    $expcorr = 1.0;
                }

                # save this all for later
                push @{$evt->{psfcorr}},  getMean($psfcorrfacts);
                push @{$evt->{psfcorrfacts}},  $psfcorrfacts;
                push @{$evt->{psfcorrtimes}},  $psfcorrtimes;
                push @{$evt->{engconv}},  $econv;
                push @{$evt->{hengconv}}, $heconv;
                push @{$evt->{sengconv}}, $seconv;
                push @{$evt->{expcorr}},  $expcorr;
                push @{$evt->{backcorr}}, $bkgexpcorr;
                push @{$evt->{exposure}}, $exposure;
                push @{$evt->{detected}}, $detected;
                push @{$evt->{orbevts}},  $orbevts;
                push @{$evt->{expomap}},  $orbexpo;
                push @{$evt->{backscl}},  $backscale;
                push @{$evt->{chisq}},    $chisq;
                push @{$evt->{dof}},      $dof;
                push @{$evt->{plindex}},  $plind;

                # remember the conversion/corrections
                $lastEconv  = $econv;
                $lastHEconv = $heconv;
                $lastSEconv = $seconv;
            }
            
            # save the source counts, bkg counts, and upper limit for
            # each band
            if ( $band eq "fb" ) {
                push @{$evt->{srccnts}}, $srccnts;
                push @{$evt->{bkgcnts}}, $bkgcnts;
            } elsif ( $band eq "hard" ) {
                push @{$evt->{hsrccnts}}, $srccnts;
                push @{$evt->{hbkgcnts}}, $bkgcnts;
            } else {
                push @{$evt->{ssrccnts}}, $srccnts;
                push @{$evt->{sbkgcnts}}, $bkgcnts;
            }

            # if this is our first time through, setup the output file
            # use the light curve as the template for output keywords
            if ( $first ) {
                $lctempl = $srclc;
                if ( !exists $fixFits{$infoFile} ) {
                    $fixFits{$infoFile} = { templ => $lctempl, gti => $allGTIFile };
                }
                $first = 0;
                ( $stat, $allfptr ) = beginLightCurveFITS( $allLCFile, $params{srcra}, $params{srcdec},
                                                           $params{trigtime}, $lctempl, 1 );
                quit( $stat ) unless $stat == 0;
            }

            # calculate the output LCs
            # at this point keep both the upper limit and the count
            # rate so we can merge upper limits 
            chat( 2, "calculating net, corrected and flux light curves\n" );
            $stat = writeUnbinnedLightCurve( $allfptr, $srclc, $bkglc, $evt, $i, $band );
            quit( $stat ) unless $stat == 0;

            # save the gti filename
            push @gtis, "$srclc\[2\]";
        }
    }
    
    # there is data in this mode if !first
    if ( !$first ) {
       
        fits_close_file( $allfptr, $stat );
        quit( $stat ) unless $stat == 0;

        # time sort
        $stat = timeSortLC( $allLCFile );
        quit( $stat ) unless $stat == 0;

        # OR the GTIs
        $stat = orGTIs( $allGTIFile, \@gtis );
        quit( $stat ) unless $stat == 0;

        # add comments unbinned LCs
        $stat = addLCComments( "$allLCFile\[RATE\]", 'src', $emin, $emax );
        quit( $stat ) unless $stat == 0;

        # save in the list to update keys
        $fixFits{$allLCFile} = { templ => $lctempl, gti => $allGTIFile };
            
        # if full band, merge phas
        if ( $band eq 'fb' ) {
            # add the background - with proper scaling
            # note that we don't bother with ratios of exposure
            # time, since they are the same
            # also note that we split here based on XRTVSUB!
            my $phanum = ( scalar( keys %arfs ) ) > 1 ? "1" : "";
            my $totalexposure = 0.0;
            foreach my $rmf ( keys %arfs ) {
                my $totbkgscl = 0.0;
                for ( my $i = 0; $i < @{$bkgexps{$rmf}}; $i++ ) {
                    $totbkgscl += $bkgexps{$rmf}[ $i ] * $bkgscls{$rmf}[ $i ];
                }
                my $avgPHA = $lcbase . "sr${phanum}.pha";
                my $bkgPHA = $lcbase . "bg${phanum}.pha";
                my $avgARF = $lcbase . "${phanum}.arf";
   
                if ( $totbkgscl > 0.0 ) {

                    my @bkgwgts = map( $totbkgexpo{$rmf} * $bkgscls{$rmf}[ $_ ] / $totbkgscl, 
                                       ( 0..$#{$bkgphas{$rmf}} ) );
                    
                    $stat = addPHAs( $bkgphas{$rmf}, \@bkgwgts, basename( $bkgPHA ),
                                     $totbkgexpo{$rmf} / $totbkgscl,
                                     'NONE', 'NONE', 'NONE' );
                    quit( $stat ) unless $stat == 0;
                    $fixFits{$bkgPHA} = { templ => $lctempl, gti => $allGTIFile };
                } else {
                    $bkgPHA = 'NONE';
                }

                # now add the source arfs
                if ( @{$arfs{$rmf}} && @{$phas{$rmf}} ) {
                    $stat = addARFs( $avgARF, $arfs{$rmf}, $exps{$rmf}, $totexpo{$rmf} );
                    quit( $stat ) unless $stat == 0;

                    $stat = addPHAs( $phas{$rmf}, undef, basename( $avgPHA ), 1.0, basename( $bkgPHA ),
                                     basename( $avgARF ), basename( $rmf ) );
                    quit( $stat ) unless $stat == 0;
                    $fixFits{$avgPHA} = { templ => $lctempl, gti => $allGTIFile };
                    $fixFits{$avgARF} = { templ => $lctempl, gti => $allGTIFile };
                    
                # don't keep the bkg pha if no source pha
                } elsif ( $bkgPHA ne 'NONE' ) {
                    push @clean, $bkgPHA;
                }
                if ( $phanum ) { $phanum++; }

                # sum the total exposure to update the combined events
                $totalexposure += $totexpo{$rmf};
            }
            # save the template for the ratio fits
            $rattempl = $lctempl;

            # if PC mode, update total events exposure
            #if ( $mode eq 'pc' ) {
            #    $stat = runSystem( genCmd( 'fparkey', { value    => $totalexposure,
            #                                            fitsfile => "$pcEvts\[EVENTS\]",
            #                                            keyword  => 'EXPOSURE' } ) );
            #    quit( $stat ) unless $stat == 0;
            #}
        }
            
        # bin them if asked to
        if ( $params{bincurve} ) {
            
            # if full band then bin, if hard band then bin
            # if soft band just apply hard band bins
            if ( $band eq 'fb' || $band eq 'hard' ) {

                # set the minimum bin time (this should be a parameter)
                my $mintime = $mode eq 'wt' ? 0.5 : 2.51;
                if ( $params{bintype} == 0 ) {
                    ( $stat, $rawBLC ) = binLCCts( $allLCFile, $params{minbincts}, 
                                                   $params{gapintv}, $mode, $mintime,
                                                   0, $params{detthresh} );
                    quit( $stat ) unless $stat == 0;
                } else {
                    ( $stat, $rawBLC ) = binLCCts( $allLCFile, $params{minbincts}, 
                                                   $params{gapintv}, $mode,
                                                   $mintime, 1, $params{detthresh},
                                                   \@{$steprates{$mode}{rates}},
                                                   \@{$steprates{$mode}{counts}} );
                    quit( $stat ) unless $stat == 0;
                }

                # cut low bins
                if ( $params{cutlowbins} || $params{"cutlowbins$mode"} ) {
                    cutLowBins( $rawBLC );
                }
            }

            # apply the bins if soft band
            if ( $band eq 'soft' ) {
                ( $stat, $rawBLC ) = applyHardToSoft( $hardBLC, $allLCFile, $mode );
            }
            
            # write the curve to disk
            chat( 2, "writing binned light curves\n" );
            $stat = writeXRTGRBLC( $rawBLC, $allBLCFile, $params{srcra}, 
                                   $params{srcdec}, $params{trigtime}, $lctempl );
            quit( $stat ) unless $stat == 0;
           
            # write the QDP file (or continue writing it)
            if ( $band eq 'fb' && $params{plotcurves} && @{$rawBLC->{time}} ) {
                genLCurveQDP( $lcplotfile, $mode, $rawBLC, $emin, $emax );
            }

            # merge the binned and unbinned into separate extensions
            $stat = runSystem( "fparkey RATEBIN \"$allBLCFile\[RATE\]\" EXTNAME" );
            quit( $stat ) unless $stat == 0;

            $stat = runSystem( "ftappend \"$allBLCFile\[RATEBIN\]\" \"$allLCFile\" chatter=1 history=no" );
            quit( $stat ) unless $stat == 0;
           
        }

        if ( $band eq 'fb' ) {
            # append the GTI 
            $stat = runSystem( "ftappend \"$allGTIFile\[1\]\" \"$allLCFile\" chatter=1 history=no" );
            quit( $stat ) unless $stat == 0;

            # add the comments posthumously
            if ( $params{bincurve} ) {
                $stat = addLCComments( "$allLCFile\[RATEBIN\]", 'src', $emin, $emax );
                quit( $stat ) unless $stat == 0;
            }
        }
           
        # remove the columns we don't want
        $stat = runSystem( genCmd( 'fdelcol', { infile  => "$allLCFile\[1\]",
                                                colname => 'BKGBACKCORR',
                                                confirm => 'no',
                                                proceed => 'yes'
                                              } ) );
        # update the checksum
        $stat = updateChecksum( $allLCFile );
        quit( $stat ) unless $stat == 0;
    }
   
    # redo for other bands, the order is important for binning, full -> hard -> soft
    # we expect there to be more counts in soft band than hard, so if we bin on
    # the hard band, there _should_ be enough counts in the soft band
    if ( $params{splitenergy} != 0.0 && !$first && $params{bincurve} ) {
        if ( $band eq 'fb' ) {
            $band = 'hard';
            $pmin = $piBrk + 1;
            $emin = $params{splitenergy};
            goto EXTRACT;
        } elsif ( $band eq 'hard' ) {
            $band = 'soft';
            $pmin = $piMin;
            $pmax = $piBrk;
            $emin = $params{minenergy};
            $emax = $params{splitenergy};

            # if hard band, save corrected curve for hardness ratio
            @{$hardBLC->{time}}       = @{$hardBLC->{timedel}} = 
            @{$hardBLC->{fracexp}}    = @{$hardBLC->{corrate}} =
            @{$hardBLC->{corrateerr}} = @{$hardBLC->{uplimit}} = 
            @{$hardBLC->{mincts}}     = ( );
            for ( my $i = 0; $i < @{$rawBLC->{time}}; $i++ ) {
                $hardBLC->{time}[ $i ]       = $rawBLC->{time}[ $i ];
                $hardBLC->{timedel}[ $i ]    = $rawBLC->{timedel}[ $i ];
                $hardBLC->{fracexp}[ $i ]    = $rawBLC->{fracexp}[ $i ];
                $hardBLC->{corrate}[ $i ]    = $rawBLC->{corrate}[ $i ];
                $hardBLC->{corrateerr}[ $i ] = $rawBLC->{corrateerr}[ $i ];
                $hardBLC->{uplimit}[ $i ]    = $rawBLC->{uplimit}[ $i ];
                $hardBLC->{mincts}[ $i ]     = $rawBLC->{mincts}[ $i ];
            }
            $hardLCFile = $allLCFile;
            goto EXTRACT;
        } elsif ( $params{bincurve} ) {

            # generate hardness ratio and save
            chat( 2, "generating hardness ratio\n" );
            $stat = genRatioFITS( $allLCFile, $hardLCFile, $ratBLCFile );

            # get the min/max for plotting
            my %stats = ( );
            ( $stat, %stats ) = fStatistic( "$ratBLCFile\[RATEBIN\]\[col HIGH=RATECOR1+ERRORCOR1\]",
                                            "HIGH" );
            if ( exists $stats{max} ) {
                $ratiomax = $stats{max} > $ratiomax ? $stats{max} : $ratiomax;
            }
            ( $stat, %stats ) = fStatistic( "$ratBLCFile\[RATEBIN\]\[col HIGH=RATECOR2+ERRORCOR2\]",
                                            "HIGH" );
            if ( exists $stats{max} ) {
                $ratiomax = $stats{max} > $ratiomax ? $stats{max} : $ratiomax;
            }
            ( $stat, %stats ) = fStatistic( "$ratBLCFile\[RATEBIN\]\[col LOW=RATECOR1-ERRORCOR1\]",
                                            "LOW" );
            if ( exists $stats{min} ) {
                $ratiomin = $stats{min} < $ratiomin ? $stats{min} : $ratiomin;
            }
            ( $stat, %stats ) = fStatistic( "$ratBLCFile\[RATEBIN\]\[col LOW=RATECOR2-ERRORCOR2\]",
                                            "LOW" );
            if ( exists $stats{min} ) {
                $ratiomin = $stats{min} < $ratiomin ? $stats{min} : $ratiomin;
            }

            # append the GTI to hard, soft, and ratio
            $stat = runSystem( "ftappend \"$allGTIFile\[1\]\" \"$allLCFile\" chatter=1 history=no" );
            quit( $stat ) unless $stat == 0;

            $stat = runSystem( "ftappend \"$allGTIFile\[1\]\" \"$hardLCFile\" chatter=1 history=no" );
            quit( $stat ) unless $stat == 0;

            $stat = runSystem( "ftappend \"$allGTIFile\[1\]\" \"$ratBLCFile\" chatter=1 history=no" );
            quit( $stat ) unless $stat == 0;

            # add comments
            $stat = addLCComments( "$allLCFile\[RATEBIN\]", 'src', $emin, $emax );
            quit( $stat ) unless $stat == 0;
            $stat = addLCComments( "$hardLCFile\[RATEBIN\]", 'src', $emin, $emax );
            quit( $stat ) unless $stat == 0;
            $stat = addLCComments( "$ratBLCFile\[RATEBIN\]", 'ratio', 
                                   $params{minenergy}, $params{maxenergy}, 
                                   $params{splitenergy} );
            quit( $stat ) unless $stat == 0;

            # save ratio for key update
            $fixFits{$ratBLCFile} = { templ => $rattempl, gti => $allGTIFile };

            # generate the qdp for plotting
            if ( $params{plotcurves} ) {
                genRatioQDP( $ratplotfile, $mode, $ratBLCFile,
                             $params{minenergy}, $params{maxenergy}, 
                             $params{splitenergy} );
            }
        } 
    
        # update energy ranges for next datamode
        $band = 'fb';
        $pmin = $piMin;
        $pmax = $piMax;
        $emin = $params{minenergy};
        $emax = $params{maxenergy};
    }

    # change mode and clear pha stuff
    if ( $mode eq 'wt' ) {
        $mode    = 'pc';
        %arfs    = ( );
        %exps    = ( );
        %phas    = ( );
        %bkgphas = ( );
        %bkgscls = ( );
        %bkgexps = ( );
        %totexpo = ( );
        %totbkgexpo = ( );
        goto EXTRACT;
    }

    # plot our qdp files
    if ( $params{plotcurves} ) {
        plotCurves( $lcplotfile );
        if ( $params{splitenergy} > 0.0 ) {
            plotRatios( $ratplotfile, $ratiomin*0.9, $ratiomax*1.1 );
        }
    }

    # write the info table
    $stat = writeInfoFile( $infoFile, \%evtData, $params{splitenergy}, \@detSources );
    quit( $stat ) unless $stat == 0;

    # fix the keywords of everything
    $stat = updateHeaders( \%fixFits );
    quit( $stat ) unless $stat == 0;

    # clean up and quit
    if ( $params{clean} ) {
        cleanup( );
    }
  
    outtro( $stat );
    return $stat;
} 
# END sub xrtgrblc


#
# getParams -
#
#       gets input parameters, expects hash ref as argument, and expects
#       that hash to have certain keys.
#
#       Inputs:
#               - hash ref of parameters
#
#       Outputs:
#               - input hash is modified
#               - returns status variable
#
sub getParams {

    my $pr     = shift;
    my $stat = 0;
    my $p;
    
    # chatter
    $p = 'chatter';
    ( $stat = getIntParam( $p, \$pr->{$p} ) ) == 0 || return $stat;

    # get the global chattiness
    libswxrtgrblc::setchat( $pr->{chatter} );

    # input files
    foreach $p ( qw( evtfiles mkffiles xhdfiles attfiles ) ) {
        ( $stat = getListParam( $p, \@{$pr->{$p}} ) ) == 0 || return $stat;
        if ( 1 == @{$pr->{$p}} ) {
            @{$pr->{$p}} = parseFitsList( @{$pr->{$p}} );
        }
        foreach my $file ( @{$pr->{$p}} ) {
            my $exists = 0;
            debug( "checking existence of $file\n" );
            $stat = fits_file_exists( $file, $exists, $stat );
            $stat = 0;
            if ( $exists == 0 ) {
                error( 1, "input file does not exist: $file\n" );
                return -1;
            }
        }
    }

    # output stuff, ra/dec, maketime expression
    foreach $p ( qw( outdir outstem object srcra srcdec psfcorrmeth maketimeexpr ) ) {
        ( $stat = getStrParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
    }

    # PC/WT source/background region lists
    foreach $p ( qw( pcreglist wtreglist ) ) {
        ( $stat = getIntParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
    }

    # Get xrtlccorr method from lccorravg parameters.
    $p = 'lccorravg';
    if(uc($pr->{psfcorrmeth}) eq 'XRTLCCORR') {
	( $stat = getBoolParam($p, \$pr->{$p} ) ) == 0 || return $stat;
    }

    # use xrtcentroid?
    $p = 'usecentroid';
    ( $stat = getBoolParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
    if ( $pr->{usecentroid} ) {
        $p = 'centroidhbox';
        ( $stat = getFltParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
    }
    
    # energies, thresholds and trigger time
    my @fltlist = qw( minenergy maxenergy splitenergy detthresh 
                      minredchisq maxredchisq pcfluxconv pclofluxconv pchifluxconv
                      wtfluxconv wtlofluxconv wthifluxconv minpccnts
                      minwtcnts wtframetime pcframetime fldsrcdist );
    foreach $p ( @fltlist ) {
        ( $stat = getFltParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
    }
    $p = 'usetrigtime';
    ( $stat = getBoolParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
    if ( $pr->{$p} ) {
        $p = 'trigtime';
        ( $stat = getFltParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
        $p = 'trigfrombat';
        ( $stat = getBoolParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
    }
   
    # misc
    foreach $p ( qw( usetargid splitorbits bincurve clean cleanspec clobber history ) ) {
        ( $stat = getBoolParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
    }
    if ( $pr->{bincurve} ) {
        $p = 'bintype';
        ( $stat = getIntParam( $p, \$pr->{$p} ) ) == 0 || return $stat;

        # counts binning
        if ( $pr->{$p} == 0 ) {
            $p = 'minbincts';
            ( $stat = getIntParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
        # step binning
        } else {
            foreach my $par ( qw( pcsteprates wtsteprates ) ) {
                $p = $par;
                ( $stat = getStrParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
            }
        }
        $p = 'gapintv';
        ( $stat = getFltParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
        $p = 'cutlowbins';
        ( $stat = getBoolParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
        $p = 'cutlowbinspc';
        ( $stat = getBoolParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
        $p = 'cutlowbinswt';
        ( $stat = getBoolParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
        $p = 'plotcurves';
        ( $stat = getBoolParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
        if ( $pr->{plotcurves} ) {
            $p = 'plotftype';
            ( $stat = getStrParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
        }
    }

    # get NH or which HI map to use
    $p = 'nhmap';
    ( $stat = getIntParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
    if ( $pr->{$p} == -1 ) {
        $p = 'usernh';
        ( $stat = getFltParam( $p, \$pr->{$p} ) ) == 0 || return $stat;
    }
    
    # do some sanity checking
    if ( scalar( @{$pr->{evtfiles}} ) != scalar( @{$pr->{mkffiles}} ) ||
         scalar( @{$pr->{evtfiles}} ) != scalar( @{$pr->{attfiles}} ) ||
         scalar( @{$pr->{evtfiles}} ) != scalar( @{$pr->{xhdfiles}} ) ) {
        error( 1, "number of event files must match number of att",
               "mkf and xhd files\n" );
        $stat = -1;
    }
    if ( !@{$pr->{evtfiles}} ) {
        error( 1, "at least one input events file must be specified\n" );
        $stat = -1;
    }
    
    if ( $pr->{outstem} eq "" || $pr->{outdir} eq "" ) {
        error( 1, "parameters outstem and outdir must be specified\n" );
        $stat = -1;
    }
    
    # check maketime expression
    if ( !$pr->{maketimeexpr} ) {
        error( 1, "maketime expression must not be empty\n" );
        $stat = -1;
    }

    # check energy ranges
    if ( $pr->{minenergy} < 0.0 || $pr->{maxenergy} <= 0.0 ) {
        error( 1, "energies must be positive definite\n" );
        $stat = -1;
    }
    if ( $pr->{minenergy} >= $pr->{maxenergy} ) {
        error( 1, "minenergy must be < maxenergy\n" );
        $stat = -1;
    }
    if ( $pr->{splitenergy} != 0.0 ) {
        if ( $pr->{splitenergy} <= $pr->{minenergy} ||
             $pr->{splitenergy} >= $pr->{maxenergy} ) {
            warnlo( 2, "parameter splitenergy outside energy range, ",
                  "assuming no split\n" );
            $pr->{splitenergy} = 0;
        }
    }
    
    my $nstat;
    
    # convert RA to degrees
    ( $nstat, $pr->{srcra} ) = convertRAStringToDegrees( $pr->{srcra} );
    $stat = $nstat unless $nstat == 0;
    
    # convert Dec to degrees
    ( $nstat, $pr->{srcdec} ) = convertDecStringToDegrees( $pr->{srcdec} );
    $stat = $nstat unless $nstat == 0;

    return $stat;
}    

#
# parseStepRates -
#
#       parses (pc|wt)steprates params. 
#       everything is global.
#       modifies %params
#
sub parseStepRates {
    
    my %steprates = (
        pc => {
            rates  => [ ],
            counts => [ ]
        },
        wt => {
            rates  => [ ],
            counts => [ ]
        }
    );
    if ( $params{bintype} == 1 ) {
        foreach my $mode ( qw( pc wt ) ) {
            foreach my $ratepair ( split /\s*,\s*/, $params{"${mode}steprates"} ) {
                my ( $rte, $cts ) = map( $_ * 1.0, split /\s*:\s*/, $ratepair );
                if ( !$rte || !$cts ) {
                    error( 1, "error parsing rate:cts pair in ",
                           "${mode}steprates: $ratepair\n" );
                    return ( -1, %steprates );
                }
                push @{$steprates{$mode}{rates}}, $rte;
                push @{$steprates{$mode}{counts}}, $cts;
            }
            my @order = sort { 
                               $steprates{$mode}{rates}[ $a ] <=> 
                               $steprates{$mode}{rates}[ $b ] 
                             } 0..$#{$steprates{$mode}{rates}};
            @{$steprates{$mode}{rates}}  = @{$steprates{$mode}{rates}}[ @order ];
            @{$steprates{$mode}{counts}} = @{$steprates{$mode}{counts}}[ @order ];
        }
        $params{minbincts} = $steprates{pc}{counts}[ 0 ] < $steprates{wt}{counts}[ 0 ] ?
                             $steprates{pc}{counts}[ 0 ] : $steprates{wt}{counts}[ 0 ];
    }
    return ( 0, %steprates );
}

#
# getCalDBFiles -
#
#       Queries the CALDB for the PC or WT mode rmfs, the teldef,
#       and the telescope alignment file
#
sub getCalDBFiles {
    
    my ( $mode, $startd, $stopd, $grade, $vsub, $telref, $alignref, $rmfref ) = @_;
    
    my $stat = 0;

    my ( $strtd, $strtt ) = split /T/, $startd;
    my ( $stpd,  $stpt  ) = split /T/, $stopd;
    
    my ( $calnam, $ext, $online, $numRet, $nFound );
    
    $calnam = $ext = $online = $numRet = $nFound = 0;
    
    my @caltypes = (
        { inst => 'XRT', type => 'TELDEF',    selec => '', ref => $telref },
        { inst => 'SC',  type => 'ALIGNMENT', selec => '', ref => $alignref }
    );
    if ( uc( $mode ) eq 'PC' ) {
        push @caltypes, { inst => 'XRT', type => 'MATRIX', 
                          selec => "datamode.eq.photon.and.grade.eq.$grade.and.xrtvsub.eq.$vsub", 
                          ref => $rmfref };
    } else {
        push @caltypes, { inst => 'XRT', type => 'MATRIX',
                          selec => "datamode.eq.windowed.and.grade.eq.$grade.and.xrtvsub.eq.$vsub", 
                          ref => $rmfref };
    }
    
    foreach my $caltype ( @caltypes ) {
        $stat = HDgtcalf( 'SWIFT', $caltype->{inst}, '-', '-', $caltype->{type}, $strtd, 
                          $strtt, $stpd, $stpt, $caltype->{selec}, 1, 1024, $calnam, 
                          $ext, $online, $numRet, $nFound, $stat );
                        
        ${$caltype->{ref}} = $calnam->[ 0 ];
        if ( $stat != 0 || $nFound == 0 ) {
            error( 1, "failed to get $caltype->{type} file from CALDB\n" );
            $stat = -1 unless $stat != 0;
            return $stat;
        }
        chat( 3, "found ${$caltype->{ref}} in CALDB\n" );

        # copy them locally - since some tasks can't handle the ftp:// etc stuff
        my $locf = catfile( $params{outdir}, basename( ${$caltype->{ref}} ) );
        if ( !-e $locf ) {
            $stat = ftcopy( ${$caltype->{ref}}, $locf, 1 );
            return $stat unless $stat == 0;
        }
        push @clean, $locf;
        ${$caltype->{ref}} = $locf;
    }

    return $stat;
}

# 
# getNH -
#
#       runs the FTool 'nh', and returns the weighted average NH.
#
sub getNH {

    my ( $ra, $dec, $eq, $map ) = @_;
    if ( !defined $eq ) {
        $eq = 2000;
    }
    if ( $map < 0 || $map > 1 ) {
        error( 1, "in getNH, map must be 0 or 1\n" );
        return ( -1, 0 );
    }
    my %nh = (
        equinox => $eq,
        ra      => $ra,
        dec     => $dec,
        size    => 3,
        disio   => 1,
        usemap  => $map,
        tchat   => 0,
        lchat   => -1
    );
    my $cmd = genCmd( 'nh', \%nh );
    if ( $map == 0 ) {
        $cmd .= " ; pget nh avwnh";
    } else {
        $cmd .= " ; pget nh alwnh";
    }
    my ( $stat, $out ) = runSystemNoChat( $cmd );
    chomp $out->[ 0 ];
    return ( $stat, $out->[ 0 ] * 1.0 );
}

#
# checkEventsFOV -
#
#       rough checks a list of event files to see if RA/Dec
#       are in FOV
#
sub checkEventsFOV {

    my ( $evtsref, $ra, $dec ) = @_;
    my $stat = 0;
    foreach my $evtfile ( @$evtsref ) {
        
        # get the X/Y of the RA/Dec
        chat( 2, "checking that RA/Dec are in FOV of $evtfile\n" );
        my ( $x, $y );
        ( $stat, $x, $y ) = skyToXY( "$evtfile+1", $ra, $dec );
        return $stat unless $stat == 0;
        
        # get the events header
        my $fptr = Astro::FITS::CFITSIO::open_file( "$evtfile+1", READONLY, $stat );
        my $head = $fptr->read_header( $stat );
        $fptr->close_file( $stat );
        return $stat unless $stat == 0;
        
        # find the bounds of the X and Y columns
        my ( $maxX, $maxY ) = ( undef, undef );
        foreach my $key ( sort keys %$head ) {
            if ( $key =~ /^TTYPE(\d+)$/ ) {
                my $num = $1;
                my $val = $head->{$key};
                $val =~ s/\s+//g;
                $val =~ s/['"]//g;
                if ( $val eq 'X' || $val eq 'x' ) {
                    $maxX = $head->{"TLMAX$num"};
                }
                if ( $val eq 'Y' || $val eq 'y' ) {
                    $maxY = $head->{"TLMAX$num"};
                }
            }
            if ( defined $maxX && defined $maxY ) { last; }
        }
        
        # do a rough check if we're within the bounds
        if ( $x < 0 || $y < 0 || $x > $maxX || $y > $maxY ) {
            error( 1, "RA/Dec outside bounds of $evtfile\n" );
            return -1;
        }
    }
    return $stat;
}


#
# createOutputDir -
#
#       creates output directory, given an input Swift event file, and a base
#       output directory. The output dir will be named OBSID_MODE, where
#       the values of OBSID and MODE are taken from the OBS_ID and DATAMODE 
#       keywords in the event file.
#       
#       Supported modes are PHOTON and WINDOWED
#
#       Inputs:
#               - fits file to read keywords from
#               - directory to start in
#               - reference to output dir
#               - reference to obsid
#               - reference to mode
#               - reference to sub-mode (e.g. w1, w2, ...)
#               - reference to middle pixel for PC mode
#               - reference to telescope PA (taken from PA_PNT keyword)
#
#       Outputs:
#               - modifies input references
#               - returns status var
#
sub createOutputDir {

    my ( $fitsFile, $dirRoot, $outDirRef, $startd, $stopd,
         $obsidref, $moderef, $submoderef, $midref, $paref, $binref,
         $graderef, $vsubref ) = @_;
   
    # get the OBS-ID and XRT mode keywords from the events file
    my $stat = 0;
    my $fptr;
    my $comment;
    
    $$moderef = 'badmode';

    chat( 4, "getting keywords from $fitsFile\n" );
    $fptr = Astro::FITS::CFITSIO::open_file( $fitsFile, READONLY, $stat );
    return $stat unless $stat == 0;
    $fptr->movabs_hdu( 1, ANY_HDU, $stat );
    $fptr->read_key_str( 'DATE-OBS', $$startd, $comment, $stat );
    $fptr->read_key_str( 'DATE-END', $$stopd, $comment, $stat );
    $fptr->read_key_str( 'OBS_ID', $$obsidref, $comment, $stat );
    $fptr->read_key_dbl( 'PA_PNT', $$paref, $comment, $stat );
    if ( $params{object} eq 'DEFAULT' || $params{object} eq 'NONE' || $params{object} eq '' ) {
        $fptr->read_key_str( 'OBJECT', $params{object}, undef, $stat );
    }

    # move to first extension
    $fptr->movabs_hdu( 2, ANY_HDU, $stat );

    # read the frame time
    $fptr->read_key_dbl( 'TIMEDEL', $$binref, $comment, $stat );
    
    $fptr->read_key_str( 'DATAMODE', $$moderef, $comment, $stat );

    # substrate voltage
    $fptr->read_key_lng( 'XRTVSUB',  $$vsubref,  undef, $stat );
    if ( $stat == KEY_NO_EXIST ) { 
        $stat = 0;
        fits_clear_errmsg( );
        warnlo( 2, "using default substrate voltage of 0 for $fitsFile\n" );
        $$vsubref = 0;
    }

    # get the grade string
    my $n = 1;
    while ( $stat != KEY_NO_EXIST ) {
        my $type;
        $fptr->read_key_str( "DSTYP$n", $type, undef, $stat );
        if ( $type =~ /GRADE/ ) {
            $fptr->read_key_str( "DSVAL$n", $$graderef, undef, $stat );
            last;
        }
    }
    fits_clear_errmsg( );
    if ( !defined $$graderef ) {
        warnlo( 2, "failed to determine grade selection for $fitsFile\n" );
    }
    
    # check the mode 
    if ( $$moderef eq 'PHOTON' ) {
        my $hlfwid;
        $$moderef = 'pc';

        # default to standard grade
        $$graderef = defined $$graderef ? "g$$graderef" : 'g0:12';

        # get the sub-mode based on field half-width
        $fptr->read_key_lng( 'WHALFWD', $hlfwid, $comment, $stat );
        if    ( $hlfwid == 245 ) { $$submoderef = 'w1'; }
        elsif ( $hlfwid == 250 ) { $$submoderef = 'w2'; }
        elsif ( $hlfwid == 300 ) { $$submoderef = 'w3'; }
        elsif ( $hlfwid == 240 ) { $$submoderef = 'w4'; }
        else {
            error( 1, "could not determine sub-mode for $fitsFile\n" );
            return -1;
        }
        $$midref = $hlfwid;
    } elsif ( $$moderef eq 'WINDOWED' ) {
        my $numcol = 0;
        $$moderef = 'wt';
        
        # default to standard grade
        $$graderef = defined $$graderef ? "g$$graderef" : 'g0:2';

        # get the sub-mode based on field half-width
        $fptr->read_key_lng( 'WMCOLNUM', $numcol, $comment, $stat );
        $$submoderef = 'w' . int( $numcol / 100 );
        if ( $$submoderef eq 'w0' ) {
            error( 1, "could not determine sub-mode for $fitsFile\n" );
            return -1;
        }
    } else {
        error( 1, "mode $$moderef not supported\n" );
        return -1;
    }
    $fptr->close_file( $stat );

    return $stat unless $stat == 0;
    
    # make out dir
    $$outDirRef = catdir( $dirRoot, $$obsidref . "_" 
                          . uc( $$moderef ) . uc( $$submoderef ) );
    if ( !-e $$outDirRef ) {
        chat( 3, "creating output dir $$outDirRef\n" );
        eval { mkpath( $$outDirRef ) };
        if ( $@ ) {
            error( 1, "failed to make output directory $$outDirRef\n" );
            $stat = -1;
        }
    } 

    return $stat;
}

#
# checkTiming -
#
#
sub checkTiming {

    my ( $evtData, $trigtime ) = @_;

    my $stat = 0;
    my $fptr;
    
    
    my @mjdrefs  = ( );
    my @timsyss  = ( );
    my @timrefs  = ( );

    # read timing keywords from each event file
    foreach my $evtf ( keys %{$evtData} ) {

        $fptr = Astro::FITS::CFITSIO::open_file( $evtf, READONLY, $stat );
        return $stat unless $stat == 0;

        $fptr->movabs_hdu( 2, ANY_HDU, $stat );
        
        my $head = $fptr->read_header( $stat );
        
        if ( exists $head->{MJDREFI} && exists $head->{MJDREFF} ) {
            # there is a precision loss here
            push @mjdrefs, $head->{MJDREFI} + $head->{MJDREFF};
        } elsif ( exists $head->{MJDREF} ) {
            push @mjdrefs, $head->{MJDREF};
        } else {
            error( 1, "MJD reference time not found in $evtf\n" );
            $stat = -1;
        }
        
        if ( exists $head->{TIMESYS} ) {
            push @timsyss, $head->{TIMESYS};
        } else {
            error( 1, "could not determine time system for $evtf\n" );
            $stat = -1;
        }

        if ( exists $head->{TIMEREF} ) {
            push @timrefs, $head->{TIMEREF};
        } else {
            error( 1, "could not determine time reference for $evtf\n" );
            $stat = -1;
        }
        
        # if no input trigger time specified, the try to get it from the event file
        if ( $params{usetrigtime} && $trigtime <= 0.0 ) {
            if ( exists $head->{TRIGTIME} ) {
                $trigtime = $head->{TRIGTIME} * 1.0;
            } 
        } elsif ( !$params{usetrigtime} && $trigtime <= 0.0 ) {
            $trigtime = $head->{TSTART};
        } elsif ( !$params{usetrigtime} && $head->{TSTART} < $trigtime ) {
            $trigtime = $head->{TSTART};
        }
    }

    # compare each one
    for ( my $i = 1; $i < @timrefs; $i++ ) {
        if ( $mjdrefs[ $i ] != $mjdrefs[ $i - 1 ] ) {
            error( 1, "mismatched MJDREFs in input event lists\n" );
            $stat = -1;
        }
        if ( $timsyss[ $i ] ne $timsyss[ $i - 1 ] ) {
            error( 1, "mismatched TIMESYS in input event lists\n" );
            $stat = -1;
        }
        if ( $timrefs[ $i ] ne $timrefs[ $i - 1 ] ) {
            error( 1, "mismatched TIMEREF in input event lists\n" );
            $stat = -1;
        }
    }

    if ( $trigtime <= 0.0 ) {
        warnlo( 2, "failed to get GRB trigger time from event lists\n" );
    }
    return ( $stat, $trigtime );
}

#
# runMakeTime -
#
#       Runs the maketime tool
sub runMakeTime {
    
    my ( $evt, $expr, $gtiRef ) = @_;
    
    my $stat = 0;
    
    $$gtiRef = catfile( $evt->{dir}, $evt->{stem} . "_orbgti.fits" );
    
    chat( 2, "running maketime to get orbits of $evt->{base}\n" );

    # make sure expr is not blank
    if ( $expr =~ /^\s*$/ ) {
        error( 1, "maketime expression empty\n",
               "failed to create GTI file for $evt->{base}\n" );
        return -1;
    }
    
    $expr =~ s/\"//g;
   
    my %maketime = (
        infile  => $evt->{mkffile},
        outfile => $$gtiRef,
        expr    => $expr,
        name    => 'anything',
        value   => 'anything',
        time    => 'TIME',
        start   => 'START',
        stop    => 'STOP',
        compact => 'no',
        copykw  => 'yes',
        histkw  => 'yes',
        prefr   => 1.0,
        postfr  => 0.0,
        clobber => 'yes'
    );
    my $cmd  = genCmd( 'maketime', \%maketime );
    $stat = runSystem( $cmd );

    push @clean, $$gtiRef;
    return $stat;
}

#
# getOrbits -
#
#       Reads the first GTI extension from a file. input list references
#       are made to hold start and stop times read from GTI.
#
sub getOrbits {

    my ( $gtiFile, $strtRef, $stopRef ) = @_;

    my $gtiBase = basename( $gtiFile );

    my ( $stat, $fptr, $nrows, $code, $width, $repeat, $col, $anynull );
    $stat = 0;

    chat( 2, "reading GTI file $gtiBase\n" );

    # open and read the START and STOP columns from the first GTI-type ext
    $fptr = Astro::FITS::CFITSIO::open_file( $gtiFile, READONLY, $stat );
    return $stat unless $stat == 0;
    $fptr->movnam_hdu( BINARY_TBL, "STDGTI", 0, $stat );
    if ( $stat != 0 ) {
        # try "GTI" instead
        $stat = 0;
        $fptr->movnam_hdu( BINARY_TBL, "GTI", 0, $stat );
        return $stat unless $stat == 0;
    }
    $fptr->get_num_rows( $nrows, $stat );
    
    $fptr->get_colnum( CASEINSEN, "START", $col, $stat );
    $fptr->get_coltype( $col, $code, $repeat, $width, $stat );
    $fptr->read_col( $code, $col, 1, 1, $nrows, 0, \@{$strtRef}, $anynull, $stat );
    
    $fptr->get_colnum( CASEINSEN, "STOP", $col, $stat );
    $fptr->get_coltype( $col, $code, $repeat, $width, $stat );
    $fptr->read_col( $code, $col, 1, 1, $nrows, 0, \@{$stopRef}, $anynull, $stat );
    
    $fptr->close_file( $stat );
    
    # time-order sort
    my @indices = sort { $strtRef->[ $a ] <=> $strtRef->[ $b ] } ( 0..$#{$strtRef} );
    @{$strtRef} = @{$strtRef}[ @indices ];
    @{$stopRef} = @{$stopRef}[ @indices ];

    return $stat;
}

#
# combinePCEventsFiles -
#
#       runs xselect to combine a list of event files
#
sub combinePCEventsFiles {

    my ( $evtRef, $gtifile, $outDir, $stem, $piMin, $piMax ) = @_;

    my $stat = 0;

    my $outf;

    chat( 2, "running xselect to combine event files\n" );

    # setup the output name 
    $outf = catfile( $outDir, $stem . "_pc_combined.evt" );
    headas_clobberfile( $outf );
    if ( -e $outf ) {
        error( 1, "failed to clobber file $outf\n" );
        return -1;
    }
    
    # get the proper dir for xselect's read function
    my $dir;
    if ( $evtRef->[ 0 ] =~ /^\// ) {
        $dir = '/';
    } else {
        $dir = './';
    }

    # write an xselect script, run it
    my $xsel = catfile( $outDir, $stem . "_pc_combine.xsel" );

    my $ret = open TMP, ">$xsel";
    if ( !$ret ) {
        error( 1, "failed to open xselect script $xsel\n" );
        return -1;
    }

    # read each events file, and combine them into one big'un
    print TMP "xsel\n";
    print TMP "set dumpcat\n";
    print TMP "read events \"$evtRef->[ 0 ]\"\n$dir\nyes\n";
    for ( my $i = 1; $i < @{$evtRef}; $i++ ) {
        print TMP "read events \"$evtRef->[ $i ]\"\n";
    }
    print TMP "filter pha_cutoff $piMin $piMax\n";
    print TMP "filter time file \"$gtifile\"\n";
    print TMP "extract events\n";
    print TMP "save events $outf\n";
    if ( -e $outf ) {
        print TMP "yes\n";
    }
    print TMP "no\n";
    print TMP "quit\nno\n";
    close TMP;
   
    # run it
    my $cmd = "xselect \@$xsel";
    $stat = runSystem( $cmd );
    
    # put script in clean list 
    push @clean, $xsel;
    
    # put combined events into cleanspec
    push @cleanspec, $outf;

    # clear xselect.log
    unlink "xselect.log";
    
    # ret the filename/status
    return ( $stat, $outf );
}

#
# detectSources -
#
sub detectSources {

    my ( $evts, $outDir, $base, $ra, $dec, $emin, $emax ) = @_;

    my $det    = catfile( $outDir, $base . "_det.fits" );
    headas_clobberfile( $det );
    if ( -e $det ) {
        error( 1, "failed to clobber file $det\n" );
        return -1;
    }
    my $tmpdet = getTmpFile( "det" );

    # run ximage to detect sources
    my $xim = catfile( $outDir, $base . "_det.xco" );
    my $ret = open TMP, ">$xim";
    if ( !$ret ) {
        error( 1, "failed to create $xim XImage script\n" );
        return -1;
    }
    print TMP "set exit_on_startfail 1;\n";
    print TMP "cey 2000;\n";
    print TMP "read size=1000 ra=$ra dec=$dec ";
    print TMP "emin=$emin emax=$emax ecol=PI \{$evts\};\n";
    print TMP "detect bright fitsdet=$det filedet=$tmpdet snr_thresh=3.5;\n";
    print TMP "exit;\n";
    close TMP;

    my $cmd = "ximage \@$xim";
    my ( $stat, $out ) = runSystem( $cmd );
    push @clean, $xim, $det;
    unlink $tmpdet;
    return ( $stat, $det );
}

#
# angdist -
#
#       calc's distance in arcsec between two ra/dec pairs
#
sub angdist {

    unless ( $#_ == 3 ) {
        return undef;
    }
    my ( $ra1, $dec1, $ra2, $dec2 ) = @_;

    my $dtor = 3.141592654 / 180.0;
    my $r0   = $ra1  * $dtor;
    my $d0   = $dec1 * $dtor;
    my $r1   = $ra2  * $dtor;
    my $d1   = $dec2 * $dtor;

    my $dst = sin( ( $d0 - $d1 ) / 2.0 )**2 + 
              cos( $d0 ) * cos( $d1 ) * sin( ( $r1 - $r0 ) / 2.0 )**2;

    if ( $dst > 1.0 ) {
        $dst = 1.0;
    } elsif ( $dst < 0.0 ) {
        $dst = 0.0;
    }
    $dst = 2.0 * asin( sqrt( $dst ) ) / $dtor;
    return $dst * 3600.0;
}

#
# getExtractionRegionsPC -
#
sub getExtractionRegionsPC {
    
    my ( $evt, $pcreglist, $ra, $dec, $pimin, $pimax, $det, $split ) = @_;

    my $start  = \@{$evt->{start}};
    my $stop   = \@{$evt->{stop}};
    my $bkgreg = \@{$evt->{bkgreg}};
    my $srcreg = \@{$evt->{srcreg}};

    # bin time for lc
    my $shtdel = 20;
    my $sblc   = "";

    # maximum encroachment (minimum safe distance)
    # on source by field source (degrees)
    my $mindist = 0.0028 * 3600.0;

    # get the SKY coords of the input position
    my ( $stat, $xpix, $ypix ) = skyToXY( "$evt->{evtsfile}+1", $ra, $dec );
    return $stat unless $stat == 0;
    
    # corresponding count rates - the higher the rate, the larger the region
    # pcreglist=1 and =2 have the same count rate bins.
    my @rates = ( 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1., 3., 5., 7., 100000. );
    
    chat( 2, "determining extraction regions for $evt->{base}\n" );

    # loop over all detected source
    my $orb  = 0;
    my $done = 0;
    foreach my $src ( @{$det} ) {
        
        chat( 3, "checking detected source at FK5(RA,Dec)=($src->{RA},$src->{DEC})\n" );
        
        # distance from source (degrees)
        my $dist = angdist( $ra, $dec, $src->{RA}, $src->{DEC} );
        chat( 3, "distance from main source is $dist arcsec.\n" );

        # check if we're to consider this a field source
        my $isGRB = 0;
        if ( $dist < $params{fldsrcdist} * 3600.0 ) {
            chat( 3, "$dist\" < fldsrcdist, assuming this is the source\n" );
            $isGRB = 1;
        } else {
            chat( 3, "$dist\" >= fldsrcdist, assuming this is a field source\n" );
            $isGRB = 0;
        }

        # if it's a field source, get the exclusion region
        if ( !$isGRB ) {
            my $fldrad = 0.0;
            my $closeness = $dist - 2.0 * $src->{HBOXSIZE};
            if ( $closeness < $mindist ) {
                $closeness = $dist - $src->{HBOXSIZE};
                if ( $closeness < $mindist ) {
                    $fldrad = $dist - $mindist;
                } else {
                    $fldrad = $src->{HBOXSIZE};
                }
            } else {
                $fldrad = 2.0 * $src->{HBOXSIZE};
            }

            # make the region string
            my $reg = sprintf( "fk5;-circle(%.5f,%.5f,%.2f\")", $src->{RA},
                               $src->{DEC}, $fldrad );
            $src->{REGION}  = $reg;
            $src->{FLDDIST} = $dist;
            chat( 3, "using exclusion region for fld src: $reg\n" );
            next;
        }

        if ( $done ) { next; }

        # if < 0.01 degrees then consider this the source in question
        for ( my $j = 0; $j < @{$start}; $j++ ) {
            $orb++;

            # extract per-orbit preliminary light curves 
            my $region = sprintf( "fk5;circle(%.5f,%.5f,40\")", $ra, $dec );
            
            $sblc = catfile( $evt->{dir}, "$evt->{stem}_orb${orb}.lc" );
            headas_clobberfile( $sblc );
            if ( -e $sblc ) {
                error( 1, "failed to clobber file $sblc\n" );
                return -1;
            }
            push @clean, $sblc;

            chat( 3, "extracting preliminary light curve from $evt->{base}\n" );
            chat( 3, "from SCC times: $start->[ $j ] -> $stop->[ $j ]\n" );

            my $cnts;
            ( $stat, $cnts ) = extractLC( $evt->{evtsfile}, "$evt->{stem}_orb${orb}", 
                                          $evt->{dir}, $pimin, $pimax, 
                                          $start->[ $j ], $stop->[ $j ], 0.0,
                                          $region, undef, $shtdel, $sblc );
            return $stat unless $stat == 0;
            
            # read the light curve so we can mess with it
            my %lc = getEmptyLC( );
            $stat = readLC( $sblc, \%lc );
            return $stat unless $stat == 0;

            # if no time bins, purge
            if ( !@{$lc{'time'}} ) {
                chat( 4, "purging $start->[$j] -> $stop->[$j], ",
                      "due to no good time bins\n" );
                splice @{$start}, $j, 1;
                splice @{$stop}, $j, 1;
                $stat = 0;
                $j--;
                next;
            }
        
            # get the max
            my $max = -1.0e20;
            my $mean = 0.0;
            for ( my $k = 0; $k < @{$lc{'time'}}; $k++ ) {
                if ( $lc{rate}[ $k ] > $max ) {
                    $max = $lc{rate}[ $k ];
                }
                $mean += $lc{rate}[ $k ];
            }
            $mean /= 1.0 * scalar( @{$lc{'time'}} );
            
            # if the max is <= 0, then force the smallest region to
            # be used, but avoid the lcstats failure
            my %lcstats;
            if ( $max <= 0.0 ) {
                $lcstats{chiprob} = 1.0;
                $lcstats{average} = 0.0;
                $lcstats{maximum} = $max;
            } else {
            
                # get the probability of constancy with lcstats
                ( $stat, %lcstats ) = lcStats( $sblc );
                # check the status - if we have one bin with low exposure
                # lcstats tanks
                if ( $stat != 0 ) {
                    $lcstats{chiprob} = 1.0;
                    $lcstats{average} = $mean;
                    $lcstats{maximum} = $max;
                }
                $stat = 0;
            }

            # get per orbit regions
            my ( $srcref, $bkgref ) = getPCRegionShapes( $evt, $pcreglist, $ra, $dec, $xpix, $ypix, $j );
            if ( !defined $srcref || !defined $bkgref ) {
                warnhi( 1, "could not determine region shapes for ",
                           "orbit $orb of $evt->{evtsfile}\n" );
                splice @{$start}, $j, 1;
                splice @{$stop}, $j, 1;
                $stat = 0;
                $j--;
                next;
            }
            my @srcRegions = @$srcref;
            my @bkgRegions = @$bkgref;
    
            # if the prob that it is constant is high (>10%), treat it as such
            my $k = $j;
            if ( $lcstats{chiprob} > 0.10 ) {
                # find the region range for the max count rate
                my $group = 0;
                while ( $group < @rates &&
                        $lcstats{average} > $rates[ $group ] ) {
                    $group++;
                }
                push @{$srcreg}, $srcRegions[ $group ];
                push @{$bkgreg}, $bkgRegions[ $group ];

            } elsif ( !$split ) {
                # find the region range for the mean count rate
                my $group = 0;
                while ( $group < @rates && 
                        $lcstats{maximum} > $rates[ $group ] ) {
                    $group++;
                }
                push @{$srcreg}, $srcRegions[ $group ];
                push @{$bkgreg}, $bkgRegions[ $group ];
            
            } else {
                for ( my $k = 0; $k < @{$lc{'time'}}; $k++ ) {
                    if ( $lc{rate}[ $k ] <= 0.0 ) {
                        splice @{$lc{'time'}}, $k, 1;
                        splice @{$lc{rate}}, $k, 1;
                        splice @{$lc{rateerr}}, $k, 1;
                        splice @{$lc{fracexp}}, $k, 1;
                        $k--;
                        next;
                    }
                }

                # slide along l.c. to get regions
                my $oldj = $j;
                my $realstart = $start->[ $j ] * 1.0;
                my $realstop  = $stop->[ $j ]  * 1.0;
                $j += slide( $start, $stop, $bkgreg, $srcreg, \%lc, 
                             $j, \@rates, \@bkgRegions, \@srcRegions );
                $j--;
                
                # fix the ends
                $start->[ $oldj ] = $realstart;
                $stop->[ $j ]     = $realstop;
            }

            # clean jitter out of the lists
            for ( my $i = $k; $i <= $j; $i++ ) {
                
                # record the actual orbit number
                push @{$evt->{orbit}}, $orb;
            
                # merge identical consecutive regions
                if ( $i > $k && $srcreg->[ $i ] eq $srcreg->[ $i - 1 ] ) {
                    chat( 3, "merging regions from $start->[ $i - 1 ] ",
                          "to $stop->[ $i ]\n" );
                    splice @{$start}, $i, 1;
                    $stop->[ $i - 1 ] = $stop->[ $i ];
                    splice @{$stop}, $i, 1;
                    splice @{$srcreg}, $i, 1;
                    splice @{$bkgreg}, $i, 1;
                    pop @{$evt->{orbit}};
                    $i--;
                    $j--;
                } 
            }
            
        }
        $done = 1;
    }
    
    for ( my $i = 0; $i < @{$start}; $i++ ) {
        my @srcstrs = split /\n/, $srcreg->[ $i ];
        my @bkgstrs = split /\n/, $bkgreg->[ $i ];
        chat( 3, "for $evt->{base} from $start->[ $i ] to $stop->[ $i ]\n" );
        chat( 3, "using source region:\n",   join "\n", @srcstrs );
        chat( 3, "and background region:\n", join "\n", @bkgstrs );
    }
    return $stat;
}

#
# getPCRegionShapes -
#
#       gets PC mode region shapes, adjusted for chip edges
#
sub getPCRegionShapes {

    my ( $evt, $pcreglist, $ra, $dec, $xpix, $ypix, $i ) = @_;

    my $stat = 0;

    # min and max allowed DET coords
    my $min = 300 - $evt->{midpix};
    my $max = 300 + $evt->{midpix};
    
    # get the source position in DET coords
    my ( $xdet, $ydet, $tmid );
   
    $tmid = ( $evt->{start}[ $i ] + $evt->{stop}[ $i ] ) / 2.0;
    
    ( $stat, $xdet, $ydet ) = pXForm( 'sky', 'det', $xpix, $ypix, $evt->{teldef}, 
                                      $evt->{align}, $evt->{attfile}, $tmid );
    return ( undef, undef ) unless $stat == 0;

    # find the distance to each chip edge
    my ( $distU, $distD, $distL, $distR );
    my @dists = ( 2.365 * ( $max - $ydet ),      # dist to top
                  2.365 * ( $ydet - $min ),      # dist to bot
                  2.365 * ( $xdet - $min ),      # dist to left
                  2.365 * ( $max - $xdet ) );    # dist to right
    my @signs = ( -1.0, 1.0, 1.0, -1.0 );       # sign of required move
    
    # sort them in increasing order
    my @sdists  = sort { $a <=> $b } @dists;
    my $mindist = $sdists[ 0 ];
    if ( $mindist < 40.0 ) {
        warnhi( 1, "source too close to chip edge in orbit "
                   . ( $i + 1 ) . " of $evt->{evtsfile}\n" );
        return ( undef, undef );
    }

    # now check the radii we want to use, depending on the pcreglist=1
    # or =2. Units are arcseconds.
    my @regions = ( );
    my @radii = ( );
    my @outii = ( );
    if($pcreglist == 2) {
	#rates:  ( 0.001, 0.005, 0.01, 0.05, 0.1, 0.5,   1.,    3.,    5.,    7., 100000.);
	@radii = (    12,    17,   22,   25,  45,  55,   65,    75,    85,    95,    105 );
	@outii = (     0,     0,    0,    0,   0,   0, 7.08, 14.16, 16.52, 18.88,  23.60 );
    } else {
	#rates:  ( 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1., 3., 5., 7., 100000. );
	@radii = (    12,    17,   22,   25,  45,  55, 65, 75, 85, 95,     105 );
	@outii = (     0,     0,    0,    0,   0,   0,  5, 10, 12, 15,      20 );
    }

    for ( my $j = 0; $j < @radii; $j++ ) {
        if ( $radii[ $j ] > $mindist ) {
            warnlo( 2, "adjusting possible $radii[ $j ]\" radius source region ",
                       "to account for chip edges\n" );
            $radii[ $j ] = $mindist;
        }
        my $regstr;
        if ( $outii[ $j ] == 0 ) {
            $regstr = sprintf( "fk5;circle(%.5f,%.5f,%d\")",
                               $ra, $dec, $radii[ $j ] );
        } else {
            $regstr = sprintf( "fk5;annulus(%.5f,%.5f,%d\",%d\")",
                               $ra, $dec, $outii[ $j ], $radii[ $j ] );
        }
        push @regions, $regstr;
    }

    # now get the corresponding background regions. Units are arcseconds.
    # Background regions are the same for both pcreglist=1 or =2.
    # move any regions so they are entirely contained within a 12' circle
    # centered at (DETX,DETY)=(300,300)
    my @bkgregions = ( );
    my @bkgradii = (  40,  40,  45,  50,  90, 110, 130, 150, 170, 190, 210 );
    my @bkgoutii = ( 140, 140, 145, 150, 190, 210, 230, 250, 270, 290, 310 );
    #my @bkgradii = (  40,  40,  40,  40,  65,  80,  90, 100, 110, 120, 130 );
    #my @bkgoutii = ( 140, 140, 140, 140, 165, 185, 190, 200, 220, 240, 260 );
    for ( my $j = 0; $j < @bkgradii; $j++ ) {

        # get the required Y-shift for the outer circle to remain in window boundaries
        my $yshift = $ydet;
        for ( my $k = 0; $k <= 1; $k++ ) {
            if ( $dists[ $k ] < $bkgoutii[ $j ] ) {
                $yshift += $signs[ $k ] * ( $bkgoutii[ $j ] - $dists[ $k ] ) / 2.365;
            }
        }
        
        # get the required X-shift for the outer circle to remain in window boundaries
        my $xshift = $xdet;
        for ( my $k = 2; $k <= 3; $k++ ) {
            if ( $dists[ $k ] < $bkgoutii[ $j ] ) {
                $xshift += $signs[ $k ] * ( $bkgoutii[ $j ] - $dists[ $k ] ) / 2.365;
            }
        }

        # now make sure it's within the 12' circle
        my $centerdist = sqrt( ( $xshift - 300 )**2.0 + ( $yshift - 300 )**2.0 ) +
                         $bkgoutii[ $j ] / 2.365;
        if ( $centerdist > 300 ) {
            my $theta = atan2( $ydet, $xdet );
            $yshift  += ( $centerdist - 300 ) * sin( $theta );
            $xshift  += ( $centerdist - 300 ) * cos( $theta );
        }

        # gen the region
        my $regstr;
        if ( $xshift == $xdet && $yshift == $ydet ) {
            $regstr = sprintf( "fk5;annulus(%.5f,%.5f,%d\",%d\")",
                               $ra, $dec, $bkgradii[ $j ], $bkgoutii[ $j ] );
        } else {
            # get the new coords
            my ( $bkgra, $bkgdec );
            ( $stat, $bkgra, $bkgdec ) = pXForm( 'det', 'sky', $xshift, $yshift, $evt->{teldef}, 
                                                 $evt->{align}, $evt->{attfile}, $tmid );
            return ( undef, undef ) unless $stat == 0;

            # write the region
            $regstr = sprintf( "fk5;circle(%.5f,%.5f,%d\")\nfk5;-circle(%.5f,%.5f,%d\")",
                               $bkgra, $bkgdec, $bkgoutii[ $j ], $ra, $dec, $bkgradii[ $j ] );
        }
        push @bkgregions, $regstr;
    }
    return ( \@regions, \@bkgregions );
}

# 
# rolldiff -
#
#       does 360 degree number line distance - always returns >= 0
#
sub rolldiff {

    my ( $ang1, $ang2 ) = @_;

    $ang1 %= 360.0;
    $ang2 %= 360.0;
    if ( $ang1 < 0.0 ) { $ang1 += 360.0; }
    if ( $ang2 < 0.0 ) { $ang2 += 360.0; }
    my $dist = abs( $ang1 - $ang2 ) % 360.0;
    if ( $dist > 180.0 ) { $dist = 360.0 - $dist; }
    return $dist;
}
    
#
# getExtractionRegionsWT -
#
sub getExtractionRegionsWT {
    
    my ( $evt, $wtreglist, $ra, $dec, $pimin, $pimax, $split ) = @_;

    my $start  = \@{$evt->{start}};
    my $stop   = \@{$evt->{stop}};
    my $bkgreg = \@{$evt->{bkgreg}};
    my $srcreg = \@{$evt->{srcreg}};

    my $stat = 0;
    
    # bin time for lc
    my $shtdel = 4;
    my $sblc   = "";

    # strip width
    my $wid = 15;

    # first get the roll angle for each orbit (it can and does change between orbits)
    my ( $fptr, $nrows, $col, $cn, $code, $rpt, $wdt, $anynull );
    my ( @timecol, @rollcol, @rollarr );
    @timecol = @rollcol = @rollarr = ( );
    
    $fptr = Astro::FITS::CFITSIO::open_file( $evt->{mkffile}, READONLY, $stat );
    return $stat unless $stat == 0;
    
    $fptr->movnam_hdu( BINARY_TBL, 'FILTER', 1, $stat );
    $fptr->get_num_rows( $nrows, $stat );
    $fptr->get_colnum( CASEINSEN, 'TIME', $cn, $stat );
    $fptr->get_coltype( $cn, $code, $rpt, $wdt, $stat );
    $fptr->read_col( $code, $cn, 1, 1, $nrows, 0, \@timecol, $anynull, $stat );
    $fptr->get_colnum( CASEINSEN, 'ROLL', $cn, $stat );
    $fptr->get_coltype( $cn, $code, $rpt, $wdt, $stat );
    $fptr->read_col( $code, $cn, 1, 1, $nrows, 0, \@rollcol, $anynull, $stat );
    $fptr->close_file( $stat );
    return $stat unless $stat == 0;
    
    # use roll in the middle! of each orbit
    my $i = 0;
    my $j = 0;
    my $tmid = ( $stop->[ $i ] + $start->[ $i ] ) / 2.0;
    foreach my $time ( @timecol ) {
        if ( $time >= $tmid ) {
            push @rollarr, $rollcol[ $j - 1 ];
            $i++;
            if ( $i >= scalar( @{$stop} ) ) {
                last;
            }
            $tmid = ( $stop->[ $i ] + $start->[ $i ] ) / 2.0;
        }
        $j++;
    }

    # adjust the orbit start/stop times so that the spacecraft is within
    # 1 degree of the roll found above
    $i = 0;
    $j = 0;
    foreach my $time ( @timecol ) {
        if ( $time >= $start->[ $i ] && 
             rolldiff( $rollcol[ $j ], $rollarr[ $i ] ) <= 1.0 ) {
            $start->[ $i ] = $timecol[ $j ];
            $i++;
            if ( $i >= scalar( @{$stop} ) ) {
                last;
            }
        }
        $j++;
    }

    # corresponding count rates - the higher the rate, the larger the
    # region.  There is a different number of count rate bins for the
    # two wtreglist=1 and =2 values.
    my @rates = ( );
    if($wtreglist == 2) {
	@rates = ( 1., 5., 10., 50., 100., 200., 300., 400., 600.0, 1000.0, 1000000. );
    } else {
	@rates = ( 1., 5., 10., 50., 100., 300., 400., 600.0, 1000.0, 1000000. );
    }

    chat( 2, "determining extraction regions for $evt->{base}\n" );

    my $orb = 0;
    for ( my $j = 0; $j < @{$start}; $j++ ) {
        
        # regions - determined empirically
        # for WT mode, the pileup can't be compensated for with an annulus,
        # so we just subtract a box around the pileup
        my $roll = $rollarr[ $orb ];
        #my @srcRegions = (
        #    sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",35\",%.1f)", $ra, $dec, $roll ),
        #    sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",71\",%.1f)", $ra, $dec, $roll ),
        #    sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",118\",%.1f)", $ra, $dec, $roll ),
        #    sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",165\",%.1f)", $ra, $dec, $roll ),
        #    sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",200\",%.1f)", $ra, $dec, $roll ),
        #    sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",236\",%.1f)\n", $ra, $dec, $roll ),
        #    sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",271\",%.1f)\n", $ra, $dec, $roll )
        #  . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",4.72\",%.1f)", $ra, $dec, $roll ),
        #    sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",306\",%.1f)\n", $ra, $dec, $roll )
        #  . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",9.43\",%.1f)", $ra, $dec, $roll ),
        #    sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",306\",%.1f)\n", $ra, $dec, $roll )
        #  . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",18.92\",%.1f)", $ra, $dec, $roll ),
        #    sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",306\",%.1f)\n", $ra, $dec, $roll )
        #  . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",37.84\",%.1f)", $ra, $dec, $roll )
        #);
	my @srcRegions = ( );
	if($wtreglist == 2) {
	    # The other case (below) was previously the only case, but this
	    # new case splits the 100-300 ct/s bin into two bins and
	    # adds a pile-up exclusion region to both.
	    @srcRegions = (
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",35\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",71\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",118\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",165\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",200\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",236\",%.1f)\n", $ra, $dec, $roll )
		. sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",4.72\",%.1f)", $ra, $dec, $roll ), #100-200 ct/s
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",260\",%.1f)\n", $ra, $dec, $roll )
		. sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",7.08\",%.1f)", $ra, $dec, $roll ), #200-300 ct/s
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",300\",%.1f)\n", $ra, $dec, $roll )
		. sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",9.43\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",450\",%.1f)\n", $ra, $dec, $roll )
		. sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",23.65\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",450\",%.1f)\n", $ra, $dec, $roll )
		. sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",28.38\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",450\",%.1f)\n", $ra, $dec, $roll )
		. sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",37.84\",%.1f)", $ra, $dec, $roll )
		);
	} else {
	    @srcRegions = (
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",35\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",71\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",118\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",165\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",200\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",236\",%.1f)\n", $ra, $dec, $roll ), #100-300 ct/s
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",300\",%.1f)\n", $ra, $dec, $roll )
		. sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",9.43\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",450\",%.1f)\n", $ra, $dec, $roll )
		. sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",23.65\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",450\",%.1f)\n", $ra, $dec, $roll )
		. sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",28.38\",%.1f)", $ra, $dec, $roll ),
		sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",450\",%.1f)\n", $ra, $dec, $roll )
		. sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",37.84\",%.1f)", $ra, $dec, $roll )
		);
	}

        # get the coords of the strip center
        my ( $xpix, $ypix, $xdet, $ydet );
       
        $tmid = ( $start->[ $j ] + $stop->[ $j ] ) / 2.0;
        
        ( $stat, $xpix, $ypix ) = skyToXY( "$evt->{evtsfile}+1", $ra, $dec );
        return $stat unless $stat == 0;
        
        ( $stat, $xdet, $ydet ) = pXForm( 'sky', 'det', $xpix, $ypix, $evt->{teldef}, 
                                          $evt->{align}, $evt->{attfile}, $tmid );
        return $stat unless $stat == 0;

        ( $stat, $xpix, $ypix ) = pXForm( 'det', 'sky', 300.5, $ydet, $evt->{teldef}, 
                                          $evt->{align}, $evt->{attfile}, $tmid );
        return $stat unless $stat == 0;


        # background regions are a little hairy for WT
        my $bkg1 = sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",450\",%.1f)\n", 
                            $xpix, $ypix, $roll );
        #my @bkgRegions = (
        #  $bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",82\",%.1f)", $ra, $dec, $roll ),
        #  $bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",118\",%.1f)", $ra, $dec, $roll ),
        #  $bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",165\",%.1f)", $ra, $dec, $roll ),
        #  $bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",212\",%.1f)", $ra, $dec, $roll ),
        #  $bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",247\",%.1f)", $ra, $dec, $roll ),
        #  $bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",283\",%.1f)", $ra, $dec, $roll ),
        #  $bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",318\",%.1f)", $ra, $dec, $roll ),
        #  $bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",353\",%.1f)", $ra, $dec, $roll ),
        #  $bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",353\",%.1f)", $ra, $dec, $roll ),
        #  $bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",353\",%.1f)", $ra, $dec, $roll ),
        #);
	my @bkgRegions = ( );
	#The exclusion regions are different for wcreglist=1 vs. =2.  
	if($wtreglist == 2) {
	    @bkgRegions = (
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",130\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",177\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",224\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",271\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",318\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",365\",%.1f)", $ra, $dec, $roll ), #100-200 ct/s
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",365\",%.1f)", $ra, $dec, $roll ), #200-300 ct/s
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",412\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",459\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",506\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",553\",%.1f)", $ra, $dec, $roll ),
		);
	} else {
	    @bkgRegions = (
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",150\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",200\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",260\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",300\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",500\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",550\",%.1f)", $ra, $dec, $roll ), #100-300 ct/s
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",640\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",1200\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",1200\",%.1f)", $ra, $dec, $roll ),
		$bkg1 . sprintf( "fk5;-rotbox(%.5f,%.5f,$wid\",1200\",%.1f)", $ra, $dec, $roll ),
	);
	}
        
        $orb++;
        
        # extract per-orbit light curves 
        my $region = sprintf( "fk5;rotbox(%.5f,%.5f,$wid\",50\",%.1f)", $ra, $dec, $roll );
        
        $sblc = catfile( $evt->{dir}, "$evt->{stem}_orb${orb}.lc" );
        headas_clobberfile( $sblc );
        if ( -e $sblc ) {
            error( 1, "failed to clobber file $sblc\n" );
            return -1;
        }
        push @clean, $sblc;

        chat( 3, "extracting preliminary light curve from $evt->{base}\n" );
        chat( 3, "from SCC times: $start->[ $j ] -> $stop->[ $j ]\n" );

        my $cnts;
        ( $stat, $cnts ) = extractLC( $evt->{evtsfile}, "$evt->{stem}_orb${orb}", 
                                      $evt->{dir}, $pimin, $pimax, 
                                      $start->[ $j ], $stop->[ $j ], 0.0,
                                      $region, undef, $shtdel, $sblc );
        return $stat unless $stat == 0;
        
        # read the light curve so we can mess with it
        my %lc = getEmptyLC( );
        $stat = readLC( $sblc, \%lc );
        return $stat unless $stat == 0;

        # if no time bins, purge
        if ( !@{$lc{'time'}} ) {
            chat( 4, "purging $start->[$j] -> $stop->[$j], ",
                  "due to no good time bins\n" );
            splice @{$start}, $j, 1;
            splice @{$stop}, $j, 1;
            $stat = 0;
            $j--;
            next;
        }
        
        # get the max
        my $max = -1.0e20;
        my $mean = 0.0;
        for ( my $k = 0; $k < @{$lc{'time'}}; $k++ ) {
            if ( $lc{rate}[ $k ] > $max ) {
                $max = $lc{rate}[ $k ];
            }
            $mean += $lc{rate}[ $k ];
        }
        $mean /= 1.0 * scalar( @{$lc{'time'}} );
        
        # if the max is <= 0, then force the smallest region to
        # be used, but avoid the lcstats failure
        my %lcstats;
        if ( $max <= 0.0 ) {
            $lcstats{chiprob} = 1.0;
            $lcstats{average} = 0.0;
            $lcstats{maximum} = $max;
        } else {
        
            # get the probability of constancy with lcstats
            ( $stat, %lcstats ) = lcStats( $sblc );
            # check the status - if we have one bin with low exposure
            # lcstats tanks
            if ( $stat != 0 ) {
                $lcstats{chiprob} = 1.0;
                $lcstats{average} = $mean;
                $lcstats{maximum} = $max;
            }
            $stat = 0;
        }

        # if the prob that it is constant is low, treat it as such
        my $k = $j;
        if ( $lcstats{chiprob} > 0.10 ) {
            # find the region range for the max count rate
            my $group = 0;
            while ( $group < @rates &&
                    $lcstats{average} > $rates[ $group ] ) {
                $group++;
            }
            push @{$srcreg}, $srcRegions[ $group ];
            push @{$bkgreg}, $bkgRegions[ $group ];

        } elsif ( !$split ) {
            # find the region range
            my $group = 0;
            while ( $group < @rates && 
                    $lcstats{maximum} > $rates[ $group ] ) {
                $group++;
            }
            push @{$srcreg}, $srcRegions[ $group ];
            push @{$bkgreg}, $bkgRegions[ $group ];
            
        # otherwise, use the slide method
        } else {
            for ( my $k = 0; $k < @{$lc{'time'}}; $k++ ) {
                if ( $lc{rate}[ $k ] <= 0.0 ) {
                    splice @{$lc{'time'}}, $k, 1;
                    splice @{$lc{rate}}, $k, 1;
                    splice @{$lc{rateerr}}, $k, 1;
                    splice @{$lc{fracexp}}, $k, 1;
                    $k--;
                    next;
                }
            }

            # slide along l.c. to get regions
            my $oldj = $j;
            my $realstart = $start->[ $j ] * 1.0;
            my $realstop  = $stop->[ $j ]  * 1.0;
            $j += slide( $start, $stop, $bkgreg, $srcreg, \%lc, 
                         $j, \@rates, \@bkgRegions, \@srcRegions );
            $j--;
            
            # fix the ends
            $start->[ $oldj ] = $realstart;
            $stop->[ $j ]     = $realstop;
        }
        
        # clean jitter out of the lists
        for ( my $i = $k; $i <= $j; $i++ ) {
            push @{$evt->{pa_pnts}}, $roll;
            push @{$evt->{orbit}}, $orb;
            if ( $i > $k && $srcreg->[ $i ] eq $srcreg->[ $i - 1 ] ) {
                chat( 3, "merging regions from $start->[ $i - 1 ] ",
                      "to $stop->[ $i ]\n" );
                splice @{$start},  $i, 1;
                $stop->[ $i - 1 ] = $stop->[ $i ];
                splice @{$stop},   $i, 1;
                splice @{$srcreg}, $i, 1;
                splice @{$bkgreg}, $i, 1;
                pop @{$evt->{pa_pnts}};
                pop @{$evt->{orbit}};
                $i--;
                $j--;
            }
        }
    }

    # print some info
    for ( my $i = 0; $i < @{$start}; $i++ ) {
        my @srcstrs = split /\n/, $srcreg->[ $i ];
        my @bkgstrs = split /\n/, $bkgreg->[ $i ];
        chat( 3, "using source region:\n",   join "\n", @srcstrs );
        chat( 3, "and background region:\n", join "\n", @bkgstrs );
        chat( 3, "for $evt->{base} from $start->[ $i ] to $stop->[ $i ]\n" );
    }
            
    return $stat;
}

#
# slide -
#
#       slides along a light curve, choosing extraction regions based on a set
#       of input rate boundaries (groups). It does this by finding the weighted 
#       average of a minimum of 20 bins (probably needs to be more robust). if 
#       the mean does not cross a rate boundary, then more bins are added until 
#       it does.
#
sub slide {
    
    my ( $start, $stop, $bgreg, $srreg, $lc, 
         $idx, $rates, $bkgRegions, $srcRegions ) = @_;
    
    # use a time-sliding cell to get extraction regions
    my $numbins  = @{$lc->{'time'}};
    my $group    = 0;
    my $maxgroup = 0;
    my $oldgroup = undef;
    my $lo       = 0;
    my $newlo    = 0;
    my $minbins  = ( $numbins < 20 ) ? $numbins : 20;
    my $hi       = ( $numbins < $minbins ) ? $numbins : $minbins;
    my $sum      = 0.0;
    my $weight   = 0.0;
    my $max      = -1.0;
    my $mean     = 0.0;
    my $splices  = 0;
    my $splclen  = 1;
    my $ntstrt   = 0.0;
    my $ntstop   = 0.0;
    
    while ( $hi <= $numbins ) {
        
        # get the weighted cell avg (weighted towards higher rates)
        for ( my $k = $newlo; $k < $hi; $k++ ) {
            my $wgt = $lc->{rateerr}[ $k ]**2.0;
            $sum    += $lc->{rate}[ $k ] * $wgt;
            $weight += $wgt;
        }
        if ( $weight <= 0.0 ) {
            $mean = 0.0;
        } else {
            $mean = $sum / $weight;
        }
        
        # find the region range
        $group = 0;
        while ( $group < @$rates && $mean > $rates->[ $group ] ) {
            $group++;
        }

        # if the mean l.c. is outside the previous range, start a new group
        if ( defined $oldgroup && $group != $oldgroup ) {

            # slide back to the earliest range violator
            # and include in new group
            my $slide = $hi - 1;
            if ( $group == 0 ) {
                while ( $slide - $lo >= $minbins &&
                        $lc->{rate}[ $slide ] <= $rates->[ $group ] ) {
                    $slide--;
                }
            } else {
                while ( $slide - $lo >= $minbins &&
                        $lc->{rate}[ $slide ] <= $rates->[ $group ] &&
                        $lc->{rate}[ $slide ] > $rates->[ $group - 1 ] ) {
                    $slide--;
                }
            }
            
            # take care of boundary cases
            my $chk = $slide + $minbins;

            # if there are less than half minbins left, just merge
            # and get the right range
            if ( $chk > $numbins ) { # && $chk - $numbins < $minbins / 2 ) {
                $max = -1.0;
                $sum = 0.0;
                $weight = 0.0;
                for ( my $k = $lo; $k < $numbins; $k++ ) {
                    my $wgt = $lc->{rateerr}[ $k ]**2.0;
                    $sum    += $lc->{rate}[ $k ] * $wgt;
                    $weight += $wgt;
                }
                if ( $weight == 0.0 ) {
                    $mean = 0.0;
                } else {
                    $mean = $sum / $weight;
                }
                $group = 0;
                while ( $group < @$rates && $mean > $rates->[ $group ] ) {
                    $group++;
                }
                $group = $group >= $oldgroup ? $group : $oldgroup;

                # splice in the new times
                $ntstrt = $splices == 0 ? $start->[ $idx ] : 
                                          $stop->[ $idx + $splices - 1 ];
                $ntstop = $lc->{'time'}[ $numbins - 1 ] + $lc->{tzero}
                          + $lc->{deltat} / 2.0;
                
                splice @{$start}, $idx + $splices, $splclen, $ntstrt;
                splice @{$stop}, $idx + $splices, $splclen, $ntstop;
                
                push @{$bgreg}, $bkgRegions->[ $group ];
                push @{$srreg}, $srcRegions->[ $group ];
                $splices++;
                last;
            
            }
     
            # splice in the new times replacing the original start/stop
            # if this is the first time through
            $ntstrt = $splices == 0 ? $start->[ $idx ] : 
                                      $stop->[ $idx + $splices - 1 ];
            $ntstop = $lc->{'time'}[ $slide ] + $lc->{tzero}
                      + $lc->{deltat} / 2.0;

            splice @{$start}, $idx + $splices, $splclen, $ntstrt;
            splice @{$stop}, $idx + $splices, $splclen, $ntstop;
            
            push @{$bgreg}, $bkgRegions->[ $oldgroup ];
            push @{$srreg}, $srcRegions->[ $oldgroup ];
            
            $splclen  = 0;
            $splices++;

            $lo       = $slide + 1;
            $hi       = $chk + 1;
            $newlo    = $lo;
            $sum      = 0.0;
            $weight   = 0.0;
            $max      = -1.0;
            $oldgroup = undef;

        # if the group is the same, then just add the next data point
        } else {
            $sum -= $lc->{rate}[ $hi - $minbins ] * $lc->{rateerr}[ $hi - $minbins ]**2.0;
            $weight -= $lc->{rateerr}[ $hi - $minbins ]**2.0;
            $newlo = $hi;
            $hi++;
            $oldgroup = $group;

            # except if we're done
            if ( $hi > $numbins ) {
                $ntstrt = $splices == 0 ? $start->[ $idx ] : 
                                          $stop->[ $idx + $splices - 1 ];
                $ntstop = $lc->{'time'}[ $numbins - 1 ] + $lc->{tzero}
                          + $lc->{deltat} / 2.0;
                
                splice @{$start}, $idx + $splices, $splclen, $ntstrt;
                splice @{$stop}, $idx + $splices, $splclen, $ntstop;
                
                push @{$bgreg}, $bkgRegions->[ $group ];
                push @{$srreg}, $srcRegions->[ $group ];
                
                $splices++;
                last;
            }
        }
    }
    return $splices;
}

#
# getSpecCorrections -
#
#       tries to fit the input spectrum ($grp) with a double absorbed power-law,
#       using both a psf corrected arf, and a non-psf corrected arf. The ratio of
#       the normalizations gives the psf correction for this spectrum. The energy 
#       conversion is found by setting the column densities of the wabs absorbing
#       columns to zero, and computing the model flux. Thus the flux returned is
#       the intrinsic flux of the source - except that redshift is not accounted 
#       for.
#
sub getSpecCorrections {
    
    my ( $dir, $stem, $grp, $bkg, $arfp, $arfnp, $rmf, $nh, $emin, $emax, $espl ) = @_;

    my $stat = 0;
    my $illformedfit = 0;
    
    my ( $xspec, @norms, $rate, $rate1, $rate2, %stats );

    # original normalization - used as a check that fit succeeded
    my $norig = 0.001;

    chat( 2, "determining spatial correction and flux conversion from $grp\n" );

    foreach my $arf ( $arfp, $arfnp ) {

        if ( $arf eq $arfnp ) {
            $xspec = catfile( $dir, "${stem}_plfit_nopsf.xcm" );
        } else {
            $xspec = catfile( $dir, "${stem}_plfit_psf.xcm" );
        }
        push @clean, $xspec;

        my $ret = open XSPC, ">$xspec";
        if ( !$ret ) {
            error( 1, "failed to open xspec file $xspec\n" );
            return -1;
        }
        print XSPC "autosave off\n";
        print XSPC "data $grp\n";
        print XSPC "response $rmf\n";
        print XSPC "arf $arf\n";
        print XSPC "backgrnd $bkg\n";
        print XSPC "ignore bad\n";
        printf XSPC "ignore \*\*-%.2f %.2f-\*\*\n", $emin, $emax;
        print XSPC "model wabs(wabs(po)) \& $nh,-1 \& 0.01 \& 2 \& $norig\n";
        print XSPC "renorm\n";
#        print XSPC "fit 10000\n";

	# Catch the ill-formed fit error and exit Xspec to avoid Xspec
	# giving an interactive prompt and making xrtgrblc appear to "hang"
	# indefinitely.

	print XSPC "if { [catch {fit 10000}] }  {\n";
	print XSPC "     puts \"MYILLFORMEDFIT = 1.  Xspec is giving up.\"\n";
	print XSPC "	 exit 1\n";
	print XSPC "}\n";

        print XSPC "tclout param 4\n";
        print XSPC "scan \$xspec_tclout \"\%f\" norm\n";
        print XSPC "puts \"MYNORM = \$norm\"\n";
        print XSPC "tclout param 3\n";
        print XSPC "scan \$xspec_tclout \"\%f\" plind\n";
        print XSPC "puts \"MYPLINDEX = \$plind\"\n";
        print XSPC "tclout dof\n";
        print XSPC "scan \$xspec_tclout \"\%f\" degof\n";
        print XSPC "puts \"MYDOF = \$degof\"\n";
        print XSPC "tclout stat\n";
        print XSPC "scan \$xspec_tclout \"\%f\" chisq\n";
        print XSPC "puts \"MYCHISQ = \$chisq\"\n";
        print XSPC "newpar 1 0\n";
        print XSPC "newpar 2 0\n";
        print XSPC "flux $emin $emax\n";
        print XSPC "tclout flux\n";
        print XSPC "scan \$xspec_tclout \"\%f\" flux\n";
        print XSPC "puts \"MYFLUX = \$flux\"\n";
        if ( $espl > 0.0 ) {
            print XSPC "flux $emin $espl\n";
            print XSPC "tclout flux\n";
            print XSPC "scan \$xspec_tclout \"\%f\" flux\n";
            print XSPC "puts \"MYFLUX1 = \$flux\"\n";
            print XSPC "flux $espl $emax\n";
            print XSPC "tclout flux\n";
            print XSPC "scan \$xspec_tclout \"\%f\" flux\n";
            print XSPC "puts \"MYFLUX2 = \$flux\"\n";
        } else {
            print XSPC "puts \"MYFLUX1 = 0.0\"\n";
            print XSPC "puts \"MYFLUX2 = 0.0\"\n";
        }

        print XSPC "tclout rate 1\n";
        print XSPC "scan \$xspec_tclout \"\%f\" rate\n";
        print XSPC "puts \"MYRATE = \$rate\"\n";
        if ( $espl > 0.0 ) {
            printf XSPC "ignore \*\*-%.2f %.2f-\*\*\n", $emin, $espl;
            print XSPC "tclout rate 1\n";
            print XSPC "scan \$xspec_tclout \"\%f\" rate1\n";
            print XSPC "if { \$rate1 > 0.0 } {\n";
            print XSPC "puts \"MYRATE1 = \$rate1\"\n";
            print XSPC "} else {\n";
            print XSPC "puts \"MYRATE1 = \$rate\"\n";
            print XSPC "}\n";
            print XSPC "notice all\n";
            print XSPC "ignore bad\n";
            printf XSPC "ignore \*\*-%.2f %.2f-\*\*\n", $espl, $emax;
            print XSPC "tclout rate 1\n";
            print XSPC "scan \$xspec_tclout \"\%f\" rate2\n";
            print XSPC "if { \$rate2 > 0.0 } {\n";
            print XSPC "puts \"MYRATE2 = \$rate2\"\n";
            print XSPC "} else {\n";
            print XSPC "puts \"MYRATE2 = \$rate\"\n";
            print XSPC "}\n";
        } else {
            print XSPC "puts \"MYRATE1 = \$rate\"\n";
            print XSPC "puts \"MYRATE2 = \$rate\"\n";
        }
        print XSPC "quit\ny\n";
        close XSPC;

        my $cmd = "xspec - $xspec";
        
        my $out;
        
        ( $stat, $out ) = runSystem( $cmd );
        return $stat unless $stat == 0;

        foreach my $line ( @{$out} ) {
            if ( $line =~ /^MYNORM =\s+([\d\.eE\+-]+)$/ ) {
                push @norms, $1 * 1.0;
            }
            # use non-psf corrected values for flux conversion
            if ( $arf eq $arfnp ) {
                if ( $line =~ /^MYPLINDEX =\s+/ ) {
                    my @splitl = split /\s+/, $line;
                    $stats{plind}   = $splitl[ 2 ] * 1.0;
                }
                if ( $line =~ /^MYFLUX =\s+([\d\.eE\+-]+)$/ ) {
                    $stats{econv} = $1 * 1.0;
                }
                if ( $line =~ /^MYFLUX1 =\s+([\d\.eE\+-]+)$/ ) {
                    $stats{seconv} = $1 * 1.0;
                }
                if ( $line =~ /^MYFLUX2 =\s+([\d\.eE\+-]+)$/ ) {
                    $stats{heconv} = $1 * 1.0;
                }
                if ( $line =~ /^MYRATE =\s+([\d\.eE\+-]+)$/ ) {
                    $rate = $1 * 1.0;
                }
                if ( $line =~ /^MYRATE1 =\s+([\d\.eE\+-]+)$/ ) {
                    $rate1 = $1 * 1.0;
                }
                if ( $line =~ /^MYRATE2 =\s+([\d\.eE\+-]+)$/ ) {
                    $rate2 = $1 * 1.0;
                }
                if ( $line =~ /^MYDOF =\s+([\d\.eE\+-]+)$/ ) {
                    $stats{dof} = $1 * 1.0;
                }
                if ( $line =~ /^MYCHISQ =\s+([\d\.eE\+-]+)$/ ) {
                    $stats{chisq} = $1 * 1.0;
                }
		if ( $line =~ /^MYILLFORMEDFIT/) {
		    $illformedfit = 1;
		}
            }
        }
    }

    # check the normalization - if renorm didn't do anything, then it failed
    if ( $norms[ 0 ] == $norig || $norms[ 1 ] == $norig ||
         $norms[ 0 ] == 0.0    || $norms[ 1 ] == 0.0 ||
	 $illformedfit == 1) {
        warnlo( 3, "power-law fit failed on $grp\n" );
        $stats{econv}   = FLT_NULL;
        $stats{heconv}  = FLT_NULL;
        $stats{seconv}  = FLT_NULL;
        $stats{psfcorrfacts} = [FLT_NULL];
        $stats{psfcorrtimes} = [FLT_NULL];
        $stats{plind}   = FLT_NULL;
        $stats{plindlo} = FLT_NULL;
        $stats{plindhi} = FLT_NULL;
        $stats{dof}     = FLT_NULL;
        $stats{chisq}   = FLT_NULL;
    } else {
        $stats{psfcorrfacts} = [$norms[ 0 ] / $norms[ 1 ]];
	$stats{psfcorrtimes} = [0.0];
        $stats{econv}  /= $rate;
        $stats{heconv} /= $rate2;
        $stats{seconv} /= $rate1;
    }
    
    return ( $stat, \%stats );
}

sub getXrtlccorrCorrection {
    my ( $mode, $srclc, $orbevts, $orbinst, $teldef, $attfile, $xhdfile, $pa_pnt, $lccorravg ) = @_;

    my $stat = 0;

    # run xrtlccorr to get psf correction
    my $outfile = getTmpFile( "fits" );
    my $outlc   = getTmpFile( "lc" );
    my $inlc    = getTmpFile( "lc" );

    my $lcstr;
    my $inststr;
    my $newinst;
    if ( uc( $mode ) eq 'PC' ) {
        $lcstr   = "$srclc\[3\]\[#row<2\]";
        $inststr = $orbinst;
    } else {
        $lcstr   = "$srclc\[3\]\[#row<3\]";
        $newinst = getTmpFile( "img" );
        $stat = ftcopy( $orbinst, $newinst, 1 );
        my $cmd = sprintf("fparkey %.12f \"$newinst\[0\]\" PA_PNT add=yes", $pa_pnt);
        $stat = runSystemNoChat( $cmd );
        $cmd = sprintf("fparkey %.12f \"$newinst\[1\]\" PA_PNT add=yes", $pa_pnt);
        $stat = runSystemNoChat( $cmd );
        $inststr = $newinst;
    }
    $stat = ftcopy( $lcstr, $inlc, 1 );
    return ($stat, [1.0], [0.0]) unless $stat == 0;
    # Set xrtlccorr parameters for assuming lccorravg=1, so no instrument map is generated.
    my $xrtlccorr = {
        lcfile         => $inlc,
        regionfile     => 'NONE',
        outfile        => $outlc,
        corrfile       => $outfile,
        teldef         => $teldef,
        aberration     => 'no',
        attinterpol    => 'no',
        attfile        => $attfile,
        srcx           => -1,
        srcy           => -1,
        psffile        => 'CALDB',
        psfflag        => 'yes',
        energy         => 1.5,
        createinstrmap => 'no',
	outinstrfile   => 'DEFAULT',
        infile         => $orbevts,
        hdfile         => $xhdfile,
        fovfile        => 'CALDB',
        vigfile        => 'CALDB',
        #nframe       => 1,
        pcnframe       => 1,
        wtnframe       => 1,
        instrfile      => $inststr,
        clobber        => 'yes',
        history        => 'no',
        chatter        => $params{chatter}
    };
    
    # If a time-dependent correction factor is desired, lccorravg=0 is set by the user,
    # and creating an instrument map is necessary.
    if($lccorravg == 0) {
	$xrtlccorr->{createinstrmap} = 'yes';
	$xrtlccorr->{wtnframe} = 10;
	$xrtlccorr->{pcnframe} = 4;
    };

    my $cmd = genCmd( 'xrtlccorr', $xrtlccorr );
    $stat = runSystem( $cmd );
    return ($stat, [1.0], [0.0]) unless $stat == 0;

    if ( defined $newinst ) {
        unlink $newinst;
    }
    unlink $inlc;

    # grab the correction factors and times
    my %corrfactstats = ( );
    my %corrtimestats = ( );
    my @times = ();
    my @corrfacts = ();

    if($lccorravg == 0) {
	# Read the TIME and CORRFACT columns of the correction factor file into arrays to use later.
	my $num_rows = 0;
	my $time_colnum = 0;
	my $corrfact_colnum = 0;
	my $time_coltype = 0;
	my $corrfact_coltype = 0;
	
	my $fp = Astro::FITS::CFITSIO::open_file($outfile, READONLY, $stat);
	$fp->movnam_hdu(ANY_HDU, 'LCCORRFACT', 0, $stat);
	$fp->get_num_rows($num_rows, $stat);
	$fp->get_colnum(CASEINSEN, 'TIME', $time_colnum, $stat);
	$fp->get_coltype($time_colnum, $time_coltype, undef, undef, $stat);
	$fp->read_col($time_coltype, $time_colnum, 1, 1, $num_rows, undef, \@times, undef, $stat);
	$fp->get_colnum(CASEINSEN, 'CORRFACT', $corrfact_colnum, $stat);
	$fp->get_coltype($corrfact_colnum, $corrfact_coltype, undef, undef, $stat);
	$fp->read_col($corrfact_coltype, $corrfact_colnum, 1, 1, $num_rows, undef, \@corrfacts, undef, $stat);
	$fp->close_file($stat);	   
    } else {
	# Read the mean TIME and mean CORRFACT into 1-element arrays.
	( $stat, %corrfactstats ) = fStatistic( $outfile, "CORRFACT" );
	( $stat, %corrtimestats ) = fStatistic( $outfile, "TIME" );
	@times = ($corrtimestats{mean});
	@corrfacts = ($corrfactstats{mean});
    }
    unlink $outfile, $outlc;

    return ( $stat, \@corrfacts, \@times );
}
    

#
# getUpperLimit -
#
#       runs Ximage's uplimit calculation on an event file for a specific region.
#       It uses the background (counts/pixel) provided, as well as the exposure map.
#       if the first argument is undefined, then it assumes you provided source and
#       background counts, and does not read any images
#
sub getUpperLimit {

    my ( $evts, $srccnts, $bkg, $sig, $base, $dir, $expo, $regf, $pmin, $pmax ) = @_;

    my $ximf;
    if ( defined $evts ) {
        $ximf = catfile( $dir, "$base.xco" );
    } else {
        $ximf = getTmpFile( "xco" );
    }
    my $ret = open TMP, ">$ximf";
    if ( !$ret ) {
        error( 1, "failed to open temporary XImage file $ximf\n" );
        return -1;
    }
    print TMP "set exit_on_startfail 1;\n";
    print TMP "cey 2000;\n";
    if ( defined $evts ) {
        print TMP "read size=800 emin=$pmin emax=$pmax ecol=PI \{$evts\}\n";
        if ( defined $expo ) {
            print TMP "read exposure size=800 \{$expo\}\n";
        }
        print TMP "uplimit background=$bkg sigma=$sig $regf\n";
    } else {
        print TMP "uplimit background=$bkg counts=$srccnts sigma=$sig\n";
    }

    # prefer the (hidden) Bayesian prior probability upper limit, but settle for 
    # classical upper limit if that method fails
    print TMP "if { \$uplimit(ullike) <= 1e-35 } {\n";
    print TMP "puts \"MYUPLIMIT = \$uplimit(value)\"\n";
    print TMP "} else {\n";
    print TMP "puts \"MYUPLIMIT = \$uplimit(ullike)\"\n";
    print TMP "}\n";
    print TMP "exit\n";
    close TMP;

    my $out;
    my $stat = 0;
    
    my $cmd = "ximage \@$ximf";

    ( $stat, $out ) = runSystem( $cmd );

    if ( defined $evts ) {
        push @clean, $ximf;
    } else {
        unlink $ximf;
    }

    return $stat unless $stat == 0;

    foreach my $line ( @{$out} ) {
        if ( $line =~ /^MYUPLIMIT =\s+([\d\.eE\+-]+)$/ ) {
            return ( $stat, $1 * 1.0 );
        }
    }
    return -1;
}
    
#
# getRegionCounts -
#
#       runs Ximage's counts task to get the total "counts" per detector pixel
#       in a region
#
sub getRegionCounts {

    my ( $img, $base, $dir, $regf ) = @_;

    my $ximf = catfile( $dir, "$base.xco" );
    my $ret = open TMP, ">$ximf";
    if ( !$ret ) {
        error( 1, "failed to open temporary XImage file $ximf\n" );
        return -1;
    }
    print TMP "set exit_on_startfail 1;\n";
    print TMP "cey 2000;\n";
    print TMP "read size=1000 \{$img\}\n";
    print TMP "counts \{$regf\}\n";
    print TMP "if { \$counts(total) <= 0.0 } {\n";
    print TMP "puts \"MYCOUNTS = 0.0\"\n";
    print TMP "} else {\n";
    print TMP "puts \"MYCOUNTS = \$counts(avgdet)\"\n";
    print TMP "}\n";
    print TMP "exit\n";
    close TMP;

    my $out;
    my $stat = 0;
    
    my $cmd = "ximage \@$ximf";

    ( $stat, $out ) = runSystem( $cmd );

    push @clean, $ximf;

    return $stat unless $stat == 0;

    my $tot;
    foreach my $line ( @{$out} ) {
        if ( $line =~ /^MYCOUNTS =\s+([\d\.eE\+-]+)$/ ) {
            $tot = $1 * 1.0;
        }
    }
    return ( $stat, $tot );
}

#
# runXRTCentroid -
#
#       runs the Swift tool xrtcentroid to refine the postion of the GRB
#
sub runXRTCentroid {

    my ( $evt, $base, $outdir, $ra, $dec, $hbox ) = @_;

    my $outf = catfile( $outdir, "${base}_xcen.txt" );

    chat( 2, "centroiding to refine coordinates\n" );

    my %xrtcentroid = ( 
        infile      => $evt,
        outfile     => "${base}_xcen.txt",
        outdir      => $outdir,
        posfile     => 'CALDB',
        calcpos     => 'yes',
        interactive => 'no',
        boxra       => $ra,
        boxdec      => $dec,
        boxradius   => $hbox,
        hist        => 'no',
        deriv       => 'no',
        clobber     => 'yes',
        chatter     => $params{chatter},
        cleanup     => 'yes'
    );

    my $cmd = genCmd( 'xrtcentroid', \%xrtcentroid );

    my $stat = runSystem( $cmd );
    return $stat unless $stat == 0;

    push @clean, $outf;

    my $ret = open FIL, "<$outf";
    if ( !$ret ) {
        error( 1, "failed to open file $outf\n" );
        return $ret;
    }
    my ( $nra, $ndec );
    while ( <FIL> ) {
        chomp;
        if ( /^RA\(degrees\) =\s+([\d\.\+-]+)/ ) {
            $nra = $1 * 1.0;
        }
        if ( /^Dec\(degrees\) =\s+([\d\.\+-]+)/ ) {
            $ndec = $1 * 1.0;
            $stat = 0;
            chat( 2, "refined GRB RA to:  ", sprintf( "%.5f\n", $nra ) );
            chat( 2, "refined GRB Dec to: ", sprintf( "%.5f\n", $ndec ) );
            last;
        }
        $stat = -1;
    }
    close FIL;
    
    return ( $stat, $nra, $ndec );
}

#
# skyToXY -
#
#       runs FTool sky2xy on file $evt, for position ($ra, $dec)
#       returns status, x-pixel and y-pixel
#
sub skyToXY {

    my ( $evt, $ra, $dec ) = @_;

    my %sky2xy = (
        infile    => $evt,
        xsky      => $ra,
        ysky      => $dec,
        xcol      => 'X',
        ycol      => 'Y',
        sensecase => 'no',
        tchat     => 10,
        lchat     => 0
    );

    my $cmd = genCmd( 'sky2xy', \%sky2xy );

    my ( $stat, $out ) = runSystem( $cmd );
    return $stat unless $stat == 0;
    foreach my $line ( @$out ) {
        chomp $line;
        if ( $line =~ /^ Output pixel coordinates:\s+(\S+?),\s*(\S+)$/ ) {
            return ( $stat, $1 * 1.0, $2 * 1.0 );
        }
    }
    return ( -1, -1, -1 );
}

#
# xyToSky -
#
#       runs FTool xy2sky on file $evt, for position ($xpix, $ypix)
#       returns status, ra and dec. Assumes DET coords
#
sub xyToSky {

    my ( $evt, $xpix, $ypix ) = @_;

    my %xy2sky = (
        infile    => $evt,
        xpix      => $xpix,
        ypix      => $ypix,
        xcol      => 'X',
        ycol      => 'Y',
        sensecase => 'no',
        tchat     => 0,
        lchat     => 0
    );

    my $cmd = genCmd( 'xy2sky', \%xy2sky );
    $cmd .= " ; pget xy2sky xsky ysky";
    
    my ( $stat, $out ) = runSystemNoChat( $cmd );
    return $stat unless $stat == 0;
    
    chomp $out->[ 0 ];
    chomp $out->[ 1 ];
    my $ra  = $out->[ 0 ] * 1.0;
    my $dec = $out->[ 1 ] * 1.0;

    return ( $stat, $ra, $dec );
}

#
# pXForm -
#
#       runs FTool pointxform
#
sub pXForm {

    my ( $from, $to, $x, $y, $teld, $algn, $att, $time ) = @_;

    $to   = uc( $to );
    $from = uc( $from );
   
    my %pointxform = (
        from       => $from, 
        to         => $to, 
        x          => $x, 
        y          => $y, 
        teldeffile => $teld, 
        alignfile  => $algn, 
        attfile    => $att, 
        time       => $time, 
        fromworld  => 'no', 
        chatter    => 0, 
        cleanup    => 'yes'
    );

    my $cmd = genCmd( 'pointxform', \%pointxform );
   
    my ( $stat, $out ) = runSystem( $cmd );
    return $stat unless $stat == 0;

    my $xto = undef;
    my $yto = undef;
    my $patt;
    if ( $to eq 'SKY' ) {
        $patt = qr/${to}[^\[]+\[\s*([\d\.-]+),\s*([\d\.-]+)/;
    } else {
        $patt = qr/${to}\s+([\d\.-]+),\s+([\d\.-]+)/;
    }
    foreach my $line ( @{$out} ) {
        if ( $line =~ /$patt/ ) {
            $xto = $1 * 1.0;
            $yto = $2 * 1.0;
        } 
    }
    if ( !defined $xto || !defined $yto ) {
        error( 1, "failed to get source ${to}X/${to}Y\n" );
        $stat = -1;
    }
    return( $stat, $xto, $yto );
}

#
# ftcopy -
#
#       runs ftcopy with clobber=yes and history=no
#
sub ftcopy {
    my ( $file1, $file2, $copyall ) = @_;
    
    $copyall = $copyall ? 'yes' : 'no';
    my %ftcopy = (
        infile  => $file1,
        outfile => $file2,
        copyall => $copyall,
        clobber => 'yes',
        chatter => 0,
        history => 'no'
    );
    my $cmd = genCmd( 'ftcopy', \%ftcopy );
    my $stat = runSystem( $cmd );
    return $stat;
}

#
# getFltFKey -
#
#       gets a float keyword from a fits file, very specific, I know.
#
sub getFltFKey {
    my ( $file, $key ) = @_;

    my ( $fptr, $stat, $comm, $val );
    $stat = 0;
    $fptr = Astro::FITS::CFITSIO::open_file( $file, READONLY, $stat );
    $fptr->read_key_dbl( $key, $val, $comm, $stat );
    $fptr->close_file( $stat );
    return ( $stat, $val );
}

#
# groupSpectrum -
#
#       tries to group a PHA spectrum to 20 ct/bin,
#       if not enough good bins, tries 10 ct/bin
#       if still not enough, copies input PHA to
#       output PHA.
#
sub groupSpectrum {

    my ( $srcspec, $grpspec, $chatter ) = @_;

    # group to 20 ct/bin
    my %grppha = (
        infile  => $srcspec,
        outfile => $grpspec,
        chatter => $params{chatter},
        comm    => 'group min 20',
        tempc   => 'exit',
        clobber => 'yes'
    );

    my $cmd = genCmd( 'grppha', \%grppha );
    chat( 2, "grouping PHA $srcspec to 20 ct/bin\n" );
    
    my $stat = runSystem( $cmd );
    quit( $stat ) unless $stat == 0;

    # check number of bins
    my %stats = ( );
    ( $stat, %stats ) = fStatistic( "$grpspec\[SPECTRUM\]\[QUALITY!=2&&GROUPING==1\]",
                                    "QUALITY" );

    # require at least 4 good bins
    if ( exists $stats{numb} && $stats{numb} >= 4 ) {
        return $stat;
    }

    # if not, try 10 ct/bin
    $grppha{comm} = 'group min 10';

    $cmd = genCmd( 'grppha', \%grppha );
    warnlo( 2, "20 ct/bin has only $stats{numb} good bins, trying 10 ct/bin\n" );
    
    $stat = runSystem( $cmd );
    quit( $stat ) unless $stat == 0;

    ( $stat, %stats ) = fStatistic( "$grpspec\[SPECTRUM\]\[QUALITY!=2&&GROUPING==1\]",
                                    "QUALITY" );

    # require at least 4 good bins
    if ( exists $stats{numb} && $stats{numb} >= 4 ) {
        return $stat;
    }

    # if still bad, just use unbinned spectrum
    warnhi( 2, "10 ct/bin has only $stats{numb} good bins, will use unbinned spectrum\n" );
    unlink $grpspec;
    ftcopy( $srcspec, $grpspec, 1 );

    return 0;
}

#
# makeArfs -
#
#       makes a psf corrected and non-psf corrected arf for input
#       spectrum and exposure map.
#
sub makeArfs {

    my ( $srcspec, $orbexpo, $specarfp, $specarfn ) = @_;
 
    if ( !$orbexpo ) {
        $orbexpo = 'NONE';
    }

    my %xrtmkarf = (
        outfile    => $specarfp,
        rmffile    => 'CALDB',
        mirfile    => 'CALDB',
        inarffile  => 'CALDB',
        expofile   => $orbexpo,
        transmfile => 'CALDB',
        psffile    => 'CALDB',
        psfflag    => 'yes',
        vigfile    => 'CALDB',
        phafile    => $srcspec,
        srcx       => -1,
        srcy       => -1,
        clobber    => 'yes',
        history    => 'no',
        chatter    => $params{chatter}
    );
    my $cmd = genCmd( 'xrtmkarf', \%xrtmkarf );
    chat( 2, "generating PSF corrected arf for $srcspec\n" );
    my $stat = runSystem( $cmd );
    return $stat unless $stat == 0;

    $xrtmkarf{outfile}  = $specarfn;
    $xrtmkarf{psfflag}  = 'no';
    $xrtmkarf{expofile} = 'NONE';

    $cmd = genCmd( 'xrtmkarf', \%xrtmkarf );
    chat( 2, "generating non-PSF corrected arf for $srcspec\n" );
    $stat = runSystem( $cmd );

    return $stat;
}

#
# fudgePSFCorr -
#
#       estimates PSF correction by taking the mean ratio of two
#       arfs - presumably these are PSF and non-PSF corrected
#
sub fudgePSFCorr {

    my ( $specarfp, $specarfn, $orbevts ) = @_;

    warnhi( 1, "Roughly estimating PSF correction for $orbevts\n" );

    my $pasted = getTmpFile( "fits" );
    
    my $ext = "\[SPECRESP\]";
    my $eng = "\[ENERG_LO>=$params{minenergy}&&ENERG_HI<=$params{maxenergy}\]";
    
    my %ftpaste = (
        infile    => "$specarfp$ext",
        pastefile => "$specarfn$ext\[col SPECRESP2=SPECRESP\]",
        outfile   => $pasted,
        copyall   => 'no',
        clobber   => 'yes',
        chatter   => 1,
        history   => 'no'
    );

    my $cmd = genCmd( 'ftpaste', \%ftpaste );
    my $stat = runSystemNoChat( $cmd );
    return ( $stat, [1.0], [0.0] ) unless $stat == 0;
   
    my $filtered = getTmpFile( "fits" );
    
    $stat = ftcopy( "$pasted$ext$eng", $filtered, 0 );
    return ( $stat, [1.0], [0.0] ) unless $stat == 0;
        
    unlink $pasted;
   
    my %stats = ( );
    ( $stat, %stats ) = fStatistic( "$filtered\[SPECRESP\]\[col DIV=SPECRESP2/SPECRESP\]",
                                    "DIV" );
    unlink $filtered;
    return ( $stat, [$stats{mean}], [0.0] );
}

sub addPHAs {

    my ( $list, $wgts, $outfile, $backscl, $backfile, $arfile, $rmfile ) = @_;

    my $stat = 0;

    # mathpha arbitrary limits
    my $mathphaMaxExprLen = 15000; # actually it's a little bigger
    my $mathphaMaxFiles   = 1000;

    # check the input list
    if ( !@$list ) {
        warnhi( 1, "no PHAs to add in addPHAs( )\n" );
        return 0;
    }

    # Kludgy!
    my $cwd = getcwd( );
    chdir $params{outdir};

    # generate mathpha expression
    my $expr = "";
    for ( my $i = 0; $i < @$list; $i++ ) {

        # calculate the expression length if we include this pha
        my $appendstr = "";
        if ( defined $wgts ) {
            $appendstr = sprintf( "%.8f*'%s'", $wgts->[ $i ], $list->[ $i ] );
        } else {
            $appendstr = sprintf( "'%s'", $list->[ $i ] );
        }
        $expr = $expr eq "" ? $appendstr : "$expr+$appendstr";

        # see if we've butted against mathpha's limits
        if ( length( $expr ) > $mathphaMaxExprLen || $i >= $mathphaMaxFiles ) {

            # we have, so do the calculation manually
            warnlo( 2, "bypassing mathpha - too many pha files to add\n" );
            $stat = bypassMathPHA( $list, $wgts, $outfile, $backscl, 
                                   $backfile, $arfile, $rmfile );
            chdir $cwd;
            return $stat;
        }
    }
    
    # setup mathpha
    my $tmpf = "$taskName.mathpha.$$.input";
    open TMP, ">$tmpf" or die "failed to open file $tmpf\n";
    print TMP "$expr\n";
    close TMP;

    my %mathpha = (
        expr      => "\@$tmpf",
        outfil    => "\!$outfile",
        units     => 'C',
        exposure  => 'CALC',
        areascal  => 1.0,
        properr   => 'no',
        errmeth   => 'GAUSS',
        auxfiles  => 'NONE',
        backfile  => $backfile,
        backscal  => $backscl,
        corrfile  => 'NONE',
        corrscal  => 'NONE',
        arfile    => $arfile,
        rmfile    => $rmfile,
        ncomments => 0,
        phaversn  => '1.2.0',
        chatter   => 9,
        divzero   => -99,
        clobber   => 'yes'
    );

    my $cmd = genCmd( 'mathpha', \%mathpha );
    $stat = runSystem( $cmd );
    unlink $tmpf;
    if ( $stat != 0 ) {
        chdir $cwd;
        return $stat;
    }

    # update the stupid keys
    $stat = runSystem( genCmd( 'fparkey', { value => $backfile,
                                            fitsfile => "$outfile\[1\]",
                                            keyword  => 'BACKFILE',
                                            add => 'yes' } ) );
    if ( $stat != 0 ) {
        chdir $cwd;
        return $stat;
    }
    $stat = runSystem( genCmd( 'fparkey', { value => $arfile,
                                            fitsfile => "$outfile\[1\]",
                                            keyword  => 'ANCRFILE',
                                            add => 'yes' } ) );
    if ( $stat != 0 ) {
        chdir $cwd;
        return $stat;
    }
    $stat = runSystem( genCmd( 'fparkey', { value => $rmfile,
                                            fitsfile => "$outfile\[1\]",
                                            keyword  => 'RESPFILE',
                                            add => 'yes' } ) );
    if ( $stat != 0 ) {
        chdir $cwd;
        return $stat;
    }
    $stat = runSystem( genCmd( 'fparkey', { value => $backscl,
                                            fitsfile => "$outfile\[1\]",
                                            keyword  => 'BACKSCAL',
                                            add => 'yes' } ) );
    chdir $cwd;
    return $stat;
}

#
# bypassMathPHA -
#
#       gets around string length/file # limits of mathpha
#       by simply creating the PHA by hand
#       
#       This assumes GAUSSIAN errors (sqrt(N))!
#
sub bypassMathPHA {

    my ( $list, $wgts, $outfile, $backscl, $backfile, $arfile, $rmfile ) = @_;

    my $stat = 0;

    my $hdfile = basename( getTmpFile( "hdr" ) );
    my $cdfile = basename( getTmpFile( "col" ) );
    my $dafile = basename( getTmpFile( "dat" ) );

    # read each input spectrum
    my $totexpo = 0.0;
    my ( $fptr, $nrows, $col, $typ, $any );
    my %channels = ( );

    for ( my $i = 0; $i < @$list; $i++ ) {
        
        # read the channel and counts columns
        my @chan = ( );
        my @cnts = ( );
        my $expo;
        $fptr = Astro::FITS::CFITSIO::open_file( $list->[ $i ], READONLY, $stat );
        $fptr->movnam_hdu( ANY_HDU, 'SPECTRUM', 0, $stat );
        if ( !defined $nrows ) {
            $fptr->get_num_rows( $nrows, $stat );
        }
        $fptr->read_key_dbl( "EXPOSURE", $expo, undef, $stat );
        $fptr->get_colnum( CASEINSEN, 'CHANNEL', $col, $stat );
        $fptr->get_coltype( $col, $typ, undef, undef, $stat );
        $fptr->read_col( $typ, $col, 1, 1, $nrows, undef, \@chan, undef, $stat );
        $fptr->get_colnum( CASEINSEN, 'COUNTS', $col, $stat );
        $fptr->get_coltype( $col, $typ, undef, undef, $stat );
        $fptr->read_col( $typ, $col, 1, 1, $nrows, undef, \@cnts, undef, $stat );
        $fptr->close_file( $stat );
        if ( $stat != 0 ) {
            return $stat;
        }

        # calculate the weighted counts in each channel
        for ( my $j = 0; $j < $nrows; $j++ ) {
            if ( !exists $channels{$chan[ $j ]} ) {
                if ( defined $wgts ) {
                    $channels{$chan[ $j ]} = $wgts->[ $i ] * $cnts[ $j ];
                } else {
                    $channels{$chan[ $j ]} = $cnts[ $j ];
                }
            } else {
                if ( defined $wgts ) {
                    $channels{$chan[ $j ]} += $wgts->[ $i ] * $cnts[ $j ];
                } else {
                    $channels{$chan[ $j ]} += $cnts[ $j ];
                }
            }
        }

        # calculate the total exposure
        $totexpo += $expo;
    }

    # fixup the total exposure as a string (no higher precision than mathpha)
    $totexpo = sprintf( "%1.6e", $totexpo );

    # open and write the data file
    open DAT, ">$dafile" or die "failed to open file $dafile\n";
    foreach my $channel ( sort { $a <=> $b } keys %channels ) {
        # round the counts to the nearest integer (ala mathpha)
        $channels{$channel} = int( $channels{$channel} + 0.5 );
        printf DAT "%d %d 0\n", $channel, $channels{$channel};
    }
    close DAT;

    # open and write the cdfile
    open CDF, ">$cdfile" or die "failed to open file $cdfile\n";
    print CDF <<EOF;
CHANNEL         J
COUNTS          J       count
QUALITY         I
EOF
    close CDF;

    # open and write the header file
    open HDR, ">$hdfile" or die "failed to open file $hdfile\n";
    print HDR <<EOF;
HDUCLASS= 'OGIP    '           / format conforms to OGIP standard
HDUCLAS1= 'SPECTRUM'           / PHA dataset (OGIP memo OGIP-92-007)
HDUVERS1= '1.2.0   '           / Obsolete - included for backwards compatibility
HDUVERS = '1.2.0   '           / Version of format (OGIP memo OGIP-92-007)
HDUCLAS2= 'DERIVED '           / WARNING This is NOT an OGIP-approved value
HDUCLAS3= 'COUNT   '           / PHA data stored as Counts (not count/s)
CREATOR = '$taskName v$taskVersion' / s/w task which wrote this dataset
TLMIN1  =                    0 / Lowest legal channel number
TLMAX1  =                 1023 / Highest legal channel number
TELESCOP= 'SWIFT   '           / mission/satellite name
INSTRUME= 'XRT     '           / instrument/detector name
FILTER  = 'NONE    '           / filter in use
EXPOSURE= $totexpo / exposure (in seconds)
AREASCAL=         1.000000E+00 / area scaling factor
BACKFILE= '$backfile'           / associated background filename
BACKSCAL= $backscl / background file scaling factor
CORRFILE= 'NONE    '           / associated correction filename
CORRSCAL=         1.000000E+00 / correction file scaling factor
RESPFILE= '$rmfile'           / associated redistrib matrix filename
ANCRFILE= '$arfile'           / associated ancillary response filename
PHAVERSN= '1992a   '           / obsolete
DETCHANS= $nrows / total number possible channels
CHANTYPE= 'PI      '           / channel type (PHA, PI etc)
POISSERR=                    T / Poissonian errors to be assumed
STAT_ERR=                    0 / no statistical error specified
SYS_ERR =                    0 / no systematic error specified
GROUPING=                    0 / no grouping of the data has been defined
EOF
    close HDR;

    # create the fits file
    $stat = ftCreate( $cdfile, $dafile, $outfile, 'SPECTRUM', $hdfile );
    unlink $cdfile, $dafile, $hdfile;
    return $stat unless $stat == 0;

    # fixup some keywords
    # 0th extension first
    my $primedit = sprintf <<EOF;
BITPIX  =                  -32 / number of bits per data pixel
CONTENT = 'PHA SPECTRUM'       / SPECTRUM xtens (at least) present
CREATOR = '$taskName v$taskVersion' / s/w task which wrote this dataset
EOF
    
    my $primedfile = dumpListToTxt( $primedit );

    my %fthedit = (
        infile     => "$outfile\[0\]",
        keyword    => "\@$primedfile",
        operation  => 'add',
        value      => 0,
        protect    => 'no',
        longstring => 'no',
        chatter    => $params{chatter},
        history    => 'no'
    );
    $stat = runSystem( genCmd( 'fthedit', \%fthedit ) );
    unlink $primedfile;
    return $stat unless $stat == 0;

    # then the 1st ext
    my $secedit = sprintf <<EOF;
-TUNIT1
-TUNIT3
EOF
    
    my $secedfile = dumpListToTxt( $secedit );

    $fthedit{infile}  = "$outfile\[1\]";
    $fthedit{keyword} = "\@$secedfile";
    $stat = runSystem( genCmd( 'fthedit', \%fthedit ) );
    unlink $secedfile;
    return $stat;
}

sub addARFs {

    my ( $outf, $arfs, $wgts, $totwgt ) = @_;

    if ( !@$arfs ) {
        warnhi( 1, "no ARFs to add in addARFs( )\n" );
        return 0;
    }

    my @dumplist = map( sprintf( "%s %.8g\n", $arfs->[$_], $wgts->[$_]/$totwgt ),
                        ( 0..$#{$arfs} ) );

    my $tmpf = dumpListToTxt( @dumplist );
    
    my %addarf = (
        list    => "\@$tmpf",
        out_ARF => "\!$outf",
        clobber => 'yes',
        chatter => 9
    );
    my $cmd  = genCmd( 'addarf', \%addarf );
    
    my $stat = runSystem( $cmd );

    unlink $tmpf;
    
    return $stat;
}

#
# writeInfoFile -
#
#       writes info about orbits, regions, detections, etc. to a fits table
#
sub writeInfoFile {
    my ( $outf, $evtData, $split, $fldregions ) = @_;
    my $stat = 0;

    headas_clobberfile( $outf );
    if ( -e $outf ) {
        error( 1, "failed to clobber file $outf\n" );
        return -1;
    }

    # column names - ugh
    my @singles  = qw( EVTSFILE MKFFILE  XHDFILE  ATTFILE  MODE  SUBMODE  OBSID );
    my @singunit =   ( '',      '',      '',      '',      '',   '',      ''    );

    my @vectors;
    my @vecunit;
    if ( $split != 0.0 ) {
        @vectors = qw( 
            START      STOP      ORBIT      SRCCNTS    BKGCNTS    HSRCCNTS
            HBKGCNTS   SSRCCNTS  SBKGCNTS 
            ENGCONV    PSFCORR   EXPCORR    BACKSCL    BACKCORR   EXPOSURE   
            DETECTED   PLINDEX   DOF        CHISQ
        );
        @vecunit = ( 
            's',       's',      '',        'counts',  'counts',  'counts',
            'counts',  'counts', 'counts', 
            'erg/cm^2/count','', '',        '',        '',        '',
            '',        '',       '',        ''
        );
    } else {
        @vectors = qw( 
            START      STOP      ORBIT      SRCCNTS    BKGCNTS        
            ENGCONV    PSFCORR   EXPCORR    BACKSCL    BACKCORR   EXPOSURE   
            DETECTED   PLINDEX   DOF        CHISQ
        );
        @vecunit = ( 
            's',       's',      '',        'counts',  'counts',
            'erg/cm^2/count','', '',        '',        '',        '',
            '',        '',       '',        ''
        );
    }

    my @regions  = qw( SRCREG    BKGREG );
    my @regunit  =   ( '',       ''     );

    # string lengths for filenames
    my $evtlen = 0;
    my $mkflen = 0;
    my $xhdlen = 0;
    my $attlen = 0;
    
    # string lengths for region columns - ugghh
    my $srcreglen = 0;
    my $bkgreglen = 0;
    
    my %coldata = ( );
    foreach my $ttyp ( @singles, @vectors, @regions ) {
        $coldata{$ttyp} = [ ];
    }
    
    my $sortfunc = sub {
        if ( $evtData->{$a}{start}[0] == $evtData->{$b}{start}[0] ) {
            return ( $evtData->{$a}{mode} cmp $evtData->{$b}{mode} );
        }
        return $evtData->{$a}{start}[0] <=> $evtData->{$b}{start}[0];
    };
    
    foreach my $key ( sort $sortfunc keys %{$evtData} ) {
                           
        my $evt = $evtData->{$key};
       
        # get max string widths
        if ( length( basename $evt->{evtsfile} ) > $evtlen ) {
            $evtlen = length( basename $evt->{evtsfile} );
        }
        if ( length( basename $evt->{mkffile} )  > $mkflen ) {
            $mkflen = length( basename $evt->{mkffile} );
        }
        if ( length( basename $evt->{xhdfile} )  > $xhdlen ) {
            $xhdlen = length( basename $evt->{xhdfile} );
        }
        if ( length( basename $evt->{attfile} )  > $attlen ) {
            $attlen = length( basename $evt->{attfile} );
        }

        # save non-region data
        foreach my $ky ( @singles ) {
            # here we push basenames, which should not affect non-filename things
            push @{$coldata{$ky}}, ( basename $evt->{lc($ky)} ) x scalar( @{$evt->{start}} );
        }

        foreach my $ky ( @vectors ) {
            push @{$coldata{$ky}}, @{$evt->{lc($ky)}};
        }

        # figure the field width and vector length for the regions
        for ( my $i = 0; $i < @{$evt->{start}}; $i++ ) {
            my $regstr = $evt->{srcreg}[ $i ];
            chomp $regstr;
            $regstr =~ s/\n/:/g;
            if ( length( $regstr )  > $srcreglen ) {
                $srcreglen = length( $regstr );
            }
            push @{$coldata{SRCREG}}, $regstr;
            $regstr = $evt->{bkgreg}[ $i ];
            chomp $regstr;
            $regstr =~ s/\n/:/g;
            if ( length( $regstr )  > $bkgreglen ) {
                $bkgreglen = length( $regstr );
            }
            push @{$coldata{BKGREG}}, $regstr;
        }
    }

    # uggggghhhhhh
    my @ttype = ( @singles,  @vectors, @regions );
    my @tunit = ( @singunit, @vecunit, @regunit );
    my @tform;
    if ( $split != 0 ) {
        @tform = ( 
            "${evtlen}A", "${mkflen}A", "${xhdlen}A", "${attlen}A", '2A', '2A',
            '11A', 'D', 'D', 'J', 'E', 'E', 'E', 'E', 'E', 'E',
            'E', 'E', 'E', 'E', 'E', 'E', 'L', 'E', 'E', 'E',
            "${srcreglen}A", "${bkgreglen}A" 
        );
    } else {
        @tform = ( 
            "${evtlen}A", "${mkflen}A", "${xhdlen}A", "${attlen}A", '2A', '2A',
            '11A', 'D', 'D', 'J', 'E', 'E', 'E', 'E', 'E', 'E', 'E', 'E',
            'L', 'E', 'E', 'E', "${srcreglen}A", "${bkgreglen}A" 
        );
    }
    
    # create the file
    my $fptr = Astro::FITS::CFITSIO::create_file( $outf, $stat );
    return $stat unless $stat == 0;
    
    # null primary array
    $fptr->write_grphdr( 1, 8, 0, [ 0 ], 0, 1, 1, $stat );
    
    my ( $datestring, $timeref );
    Astro::FITS::CFITSIO::fits_get_system_time( $datestring, $timeref, $stat );
    $fptr->write_key_str( 'DATE', $datestring, 'File creation date', $stat );

    # add the table
    my $nrows = scalar( @{$coldata{EVTSFILE}} );
    my $nfld  = scalar( @ttype );
    $fptr->insert_btbl( $nrows, $nfld, \@ttype, \@tform, \@tunit, "XRTINFO", 0, $stat );

    my $col;
        
    my $i = 0;
    foreach my $colname ( @singles, @vectors, @regions ) {
        
        # add the columns
        $fptr->get_colnum( CASEINSEN, "$colname", $col, $stat );
        if ( $tform[ $i ] =~ /A/ ) {
            $fptr->write_col_str( $col, 1, 1, $nrows, \@{$coldata{$colname}}, $stat );
        } elsif ( $tform[ $i ] =~ /D/ ) {
            $fptr->write_col_dbl( $col, 1, 1, $nrows, \@{$coldata{$colname}}, $stat );
        } elsif ( $tform[ $i ] =~ /E/ ) {
            if ( $colname =~ /uplim/i || $colname =~ /plind/i ||
                 $colname =~ /dof/i   || $colname =~ /chisq/i ) {
                $fptr->write_colnull_flt( $col, 1, 1, $nrows, 
                                          \@{$coldata{$colname}}, FLT_NULL, $stat );
            } else {
                $fptr->write_col_flt( $col, 1, 1, $nrows, \@{$coldata{$colname}}, $stat );
            }
        } elsif ( $tform[ $i ] =~ /J/ ) {
            $fptr->write_col_lng( $col, 1, 1, $nrows, \@{$coldata{$colname}}, $stat );
        } elsif ( $tform[ $i ] =~ /L/ ) {
            $fptr->write_col_log( $col, 1, 1, $nrows, \@{$coldata{$colname}}, $stat );
        }
        $i++;
    }

    # add a field region extension with all info from Ximage detect
    @ttype = qw( SRCRATE SRCRATE_ERR X Y VIGNET RA DEC ERRRAD HBOXSIZE PROB SNR REGION );
    @tform = qw(   1E       1E      1E 1E   1E  1D  1D   1E       1E    1E   1E   80A );
    @tunit = ( 'counts/s', 'counts/s', 'pixel', 'pixel', '', 
               'deg', 'deg', 'arcsec', 'arcsec', '', '', '' );

    # check that we have all the data foreach row
    # we might have a row for the GRB, which is NOT a field source
    # and has no single associated region
    $nrows = @{$fldregions};
    for ( my $i = 0; $i < @{$fldregions}; $i++ ) {
        foreach my $colname ( @ttype ) {
            if ( !exists $fldregions->[ $i ]->{$colname} ) {
                $nrows--;
                splice @{$fldregions}, $i, 1;
            }
        }
    }
    
    # add the table, even if it will be empty
    $nfld = scalar( @ttype );
    $fptr->insert_btbl( $nrows, $nfld, \@ttype, \@tform, \@tunit, "XRTFLDREG", 0, $stat );
    
    if ( $nrows > 0 ) {
        my $i = 0;

        foreach my $colname ( @ttype ) {
            
            # add the columns
            my @data = map ( $_->{$colname}, @{$fldregions} );
            
            $fptr->get_colnum( CASEINSEN, "$colname", $col, $stat );
            
            if ( $tform[ $i ] =~ /A/ ) {
                $fptr->write_col_str( $col, 1, 1, $nrows, \@data, $stat );
            } elsif ( $tform[ $i ] =~ /D/ ) {
                $fptr->write_col_dbl( $col, 1, 1, $nrows, \@data, $stat );
            } elsif ( $tform[ $i ] =~ /E/ ) {
                $fptr->write_col_flt( $col, 1, 1, $nrows, \@data, $stat );
            }
            $i++;
        }
    } 
    $fptr->close_file( $stat );
    return $stat;

    # meh
}

#
# timeSortLC -
#
#       time sorts the light curves
#
sub timeSortLC {

    my $lcfile = shift;

    chat( 2, "time sorting unbinned light curves\n" );

    # time sort the lc using fmemsort
    my $tmpf = getTmpFile( "fits" );
    my $fmemsort = {
        infile    => "$lcfile\[RATE\]",
        outfile   => $tmpf,
        columns   => 'TIME',
        method    => 'heap',
        ascend    => 'yes',
        load2mem  => 'no',
        copyprime => 'yes',
        copyall   => 'yes',
        unique    => 'no',
        history   => 'no',
        clobber   => 'yes'
    };
    my $stat = runSystem( genCmd( 'fmemsort', $fmemsort ) );
    if ( !-e $tmpf ) {
        warnhi( 1, "failed to sort file $lcfile\[1\], check number of rows\n" );
        return 0;
    }
    if ( !move( $tmpf, $lcfile ) ) {
        error( 1, "failed to sort $lcfile\n" );
        return -1;
    }
    return 0;
}

#
# orGTIs -
#
#       runs mgtime to or a bunch of GTI files
#
sub orGTIs {
    my ( $outfile, $listref ) = @_;

    my $stat = 0;

    my @totstart = ( );
    my @totstop  = ( );

    # read the starts/stops from each file
    foreach my $file ( @$listref ) {
        my @start = ( );
        my @stop  = ( );
        $stat = getOrbits( $file, \@start, \@stop );
        return $stat unless $stat == 0;

        push @totstart, @start;
        push @totstop,  @stop;
    }

    # time sort based on start
    my @order = sort { $totstart[ $a ] <=> $totstart[ $b ] } 0..$#totstart;
    @totstart = @totstart[ @order ];
    @totstop  = @totstop[ @order ];

    # or them
    for ( my $i = 1; $i < @totstart; $i++ ) {
        if ( $totstop[ $i - 1 ] >= $totstart[ $i ] ) {
            $totstop[ $i - 1 ] = $totstop[ $i ];
            splice @totstart, $i, 1;
            splice @totstop,  $i, 1;
            $i--;
        }
    }

    # dump them to the output file
    $stat = writeSwiftMETGTI( \@totstart, \@totstop, $outfile );

    return $stat;
}
#    my @list = @$listref;
#
#    my $stat = 0;
#
#    # if we only have one, duplicate it 
#    # This is for keyword consistency
#    my $copygti;
#    if ( $#list == 0 ) {
#        $copygti = getTmpFile( "gti" );
#        $stat = ftcopy( $list[ 0 ], $copygti, 0 );
#        return $stat unless $stat == 0;
#        push @list, $copygti;
#    }
#
#    my $mgtime = {
#        merge    => 'OR',
#        instarts => 'START',
#        instops  => 'STOP',
#        indates  => 'MJDREF',
#        outstart => 'START',
#        outstop  => 'STOP'
#    };
#
#    # do only 100 at a time (well below mgtime's limit)
#    # this puts a limit of ~4000 GTIs that can be merged with this routine
#    # hopefully we'll never get there
#    my $strt = 0;
#    my $stop = 99 > $#list ? $#list : 99;
#    my @tmpgtis = ( );
#    while( $stop <= $#list && $strt < $stop ) {
#
#
#        my $tmpf = dumpListToTxt( @list[ $strt..$stop ] );
#        my $tmpout = getTmpFile( "gti" );
#        $mgtime->{ingtis} = "\@$tmpf";
#        $mgtime->{outgti} = $tmpout;
#        $stat = runSystem( genCmd( 'mgtime', $mgtime ) );
#
#        unlink $tmpf;
#       
#        push @tmpgtis, $tmpout;
#        $strt = $stop + 1;
#        $stop = ( $strt + 99 ) > $#list ? $#list : $strt + 99;
#    }
#    if ( $#tmpgtis == 0 ) {
#        if ( !move( $tmpgtis[ 0 ], $outfile ) ) {
#            error( 1, "failed to move $tmpgtis[ 0 ] to $outfile\n" );
#            return -1;
#        }
#    } else {
#        my $tmpf = dumpListToTxt( @tmpgtis );
#        $mgtime->{ingtis} = "\@$tmpf";
#        $mgtime->{outgti} = $outfile;
#        $stat = runSystem( genCmd( 'mgtime', $mgtime ) );
#
#        unlink $tmpf, @tmpgtis;
#    }
#    if ( defined $copygti ) { unlink $copygti; }
#    return $stat unless $stat == 0;
#
#    # strip the history keys
#    my $fptr = Astro::FITS::CFITSIO::open_file( $outfile, READWRITE, $stat );
#    $fptr->movnam_hdu( BINARY_TBL, "STDGTI", 0, $stat );
#    while ( $stat == 0 ) {
#        $fptr->delete_key( "HISTORY", $stat );
#    }
#    if ( $stat == KEY_NO_EXIST ) {
#        fits_clear_errmsg( );
#        $stat = 0;
#    }
#    $fptr->close_file( $stat );
#
#    return $stat;
#}

#
# writeSwiftMETGTI -
#
#       writes a GTI file in Swift MET time
#
sub writeSwiftMETGTI {
    my ( $start, $stop, $gti ) = @_;

    my $tstart = sprintf( "%.14e", $start->[ 0 ] );
    my $tstop  = sprintf( "%.14e", $stop->[ $#$stop ] );

    my $hdfile = getTmpFile( "txt" );
    my $cdfile = getTmpFile( "txt" );
    my $dtfile = getTmpFile( "txt" );

    open HDR, ">$hdfile" or die "failed to open file $hdfile\n";
    print HDR <<EOF;
HDUCLASS= 'OGIP    '            / Conforms to OGIP/GSFC standards
HDUCLAS1= 'GTI     '            / Contains good time intervals
HDUCLAS2= 'STANDARD'            / Contains standard good time intervals
HDUVERS = '1.0.0   '            / Version of GTI header
TIMEZERO=                   0.  / Zero-point offset for TIME column
MJDREF  = 5.191000074287037E+04 / MJD Epoch of TIME = 0
TSTART  = $tstart  / Start time of GTI
TSTOP   = $tstop  / Stop  time of GTI
EOF
    close HDR;

    open CDF, ">$cdfile" or die "failed to open file $cdfile\n";
    print CDF <<EOF;
START           D   s
STOP            D   s
EOF
    close CDF;

    open DAT, ">$dtfile" or die "failed to open file $dtfile\n";
    for ( my $i = 0; $i < @$start; $i++ ) {
        printf DAT "%.15e %.15e\n", $start->[ $i ], $stop->[ $i ];
    }
    close DAT;

    my $stat = ftCreate( $cdfile, $dtfile, $gti, 'STDGTI', $hdfile );
    push @clean, $hdfile, $cdfile, $dtfile;
    return $stat;
}

#
# binLCCts -
#
#       bins a light curve based on counts
#
sub binLCCts {
    my ( $lcfile, $cts, $gapintv, $mode,
         $mintime, $step, $detthr, $rates, $counts ) = @_;

    # open the unbinned light curve if we are binning a file
    # otherwise we are being called recursively 
    my $null = FLT_NULL;
    my $any = 0;
    my $stat = 0;
    my $nrows = 0;
    my $optrows = 0;
    my $fptr;
    my $recursing = 0;
    if ( ref( $lcfile ) =~ /HASH/ ) {
        $recursing = 1;
    }
    if ( !$recursing ) {
        $fptr = Astro::FITS::CFITSIO::open_file( "$lcfile\[1\]", READONLY, $stat );
        $fptr->get_num_rows( $nrows, $stat );
        $fptr->get_rowsize( $optrows, $stat );
        if ( $nrows < 1 ) { return 0; }
        if ( $optrows > $nrows - 1 ) { $optrows = $nrows - 1; }
    } else {
        $nrows = scalar( @{$lcfile->{time}} );
        if ( $nrows < 1 ) { return ( 0, $lcfile ); }
    }

    # in the unusual case that the first data point has a gap afterward,
    # we need to set the default 
    if ( $recursing ) {
        chat( 2, "trying to merge low count bins\n" );
    } elsif ( $step ) { 
        $cts = $counts->[ 0 ];
        chat( 2, "step binning light curve with min S/N $detthr\n" );
    } else {    
        chat( 2, "binning light curve to $cts counts/bin\n" );
    }

    # clone scalar light curve params, and GTI
    my %outLC = getEmptyLC( );

    # read the first chunk
    my $row = 1;
    my $i = 1;
    my %columns = ( );
    my @colnames = qw( time timedel fracexp ratenet errornet ratecor errorcor
                       rateflux errorflux rateraw errorraw ratebkg errorbkg bkgbackcorr );
    my %colmap = (
        time       => 'time',
        timedel    => 'timedel',
        fracexp    => 'fracexp',
        ratenet    => 'netrate',
        errornet   => 'netrateerr',
        ratecor    => 'corrate',
        errorcor   => 'corrateerr',
        rateflux   => 'ergrate',
        errorflux  => 'ergrateerr',
        rateraw    => 'rate',
        errorraw   => 'rateerr',
        ratebkg    => 'bkgrate',
        errorbkg   => 'bkgrateerr',
        bkgbackcorr => 'bkgbackcorr'
    );
    if ( !$recursing ) {
        debug( "reading $optrows rows from file $lcfile\[1\]\n" );
    }
    foreach my $col ( @colnames ) {
        if ( !$recursing ) {
            my $num = 0;
            $fptr->get_colnum( CASEINSEN, $col, $num, $stat );
            $columns{$col}{colnum} = $num;
            my $type;
            $fptr->get_coltype( $columns{$col}{colnum}, $type, undef, undef, $stat );
            $columns{$col}{type} = $type;
            @{$columns{$col}{val}} = ( 0.0 ) x $optrows;
            $fptr->read_col( $columns{$col}{type}, $columns{$col}{colnum}, $row, 1, 
                             $optrows + 1, $null, \@{$columns{$col}{val}}, $any, $stat );
        } else {
            $columns{$col}{val} = $lcfile->{$colmap{$col}};
        }
    }
    $row++;
    
    # bin totals
    my $timeweight = $columns{timedel}{val}[ 0 ] * $columns{fracexp}{val}[ 0 ];
    my $totcts  = $timeweight * $columns{rateraw}{val}[ 0 ];
    my $toterr  = ( $timeweight * $columns{errorraw}{val}[ 0 ] )**2.0;
    my $totncts = $timeweight * $columns{ratenet}{val}[ 0 ];
    my $totnerr = ( $timeweight * $columns{errornet}{val}[ 0 ] )**2.0;
    my $totccts = $timeweight * $columns{ratecor}{val}[ 0 ];
    my $totcerr = ( $timeweight * $columns{errorcor}{val}[ 0 ] )**2.0;
    my $totects = $timeweight * $columns{rateflux}{val}[ 0 ];
    my $toteerr = ( $timeweight * $columns{errorflux}{val}[ 0 ] )**2.0;
    my $totbcts = $timeweight * $columns{ratebkg}{val}[ 0 ];
    my $totberr = ( $timeweight * $columns{errorbkg}{val}[ 0 ] )**2.0;
    my $totbcor = $timeweight * $columns{bkgbackcorr}{val}[ 0 ];
    my $totccor = 0.0; # weighted average correction
    my $totecon = 0.0; # weighted average flux conversion
    my $totccwt = 0.0; # total weight for above
    my $upper = FLT_NULL;

    # S/N = S - scl*B / sqrt( S + scl^2 * B )
    my $bgvari   = $totbcts * $columns{bkgbackcorr}{val}[ 0 ];
    my $snr      = $totcts - $totbcts;
    if ( $snr + $bgvari > 0.0 ) {
        $snr /= sqrt( $snr + $bgvari );
    } else {
        $snr = 0.0;
    }
    
    # last bin totals
    my $oldcts  = 0.0;
    my $olderr  = 0.0;
    my $oldncts = 0.0;
    my $oldnerr = 0.0;
    my $oldccts = 0.0;
    my $oldcerr = 0.0;
    my $oldects = 0.0;
    my $oldeerr = 0.0;
    my $oldbcts = 0.0;
    my $oldberr = 0.0;
    my $oldbcor = 0.0;
    my $oldweight  = 0.0;

    # starting conditions
    my $start  = $columns{time}{val}[ 0 ] - $columns{timedel}{val}[ 0 ] / 2.0;
    my $stop   = $columns{time}{val}[ 0 ] + $columns{timedel}{val}[ 0 ] / 2.0;
    my $hasBin = 0;
    my $hasGap = 0;

    #
    # anonymous subs (for (relative) cleanliness)
    #

    #
    # add a bin
    #
    my $addbin = sub {
        # set the stop time
        if ( $i - 1 < @{$columns{time}{val}} ) {
            $stop = $columns{time}{val}[ $i - 1 ] + $columns{timedel}{val}[ $i - 1 ] / 2.0;
        } else {
            $stop = $columns{time}{val}[ $i - 2 ] + $columns{timedel}{val}[ $i - 2 ] / 2.0;
        }

        # push the time, rate, etc.
        push @{$outLC{timedel}},     $stop - $start;
        push @{$outLC{time}},        ( $start + $stop ) / 2.0;
        push @{$outLC{rate}},        $totcts / $timeweight;
        push @{$outLC{rateerr}},     sqrt( $toterr ) / $timeweight;
        push @{$outLC{netrate}},     $totncts / $timeweight;
        push @{$outLC{netrateerr}},  sqrt( $totnerr ) / $timeweight;
        push @{$outLC{corrate}},     $totccts / $timeweight;
        push @{$outLC{corrateerr}},  sqrt( $totcerr ) / $timeweight;
        push @{$outLC{ergrate}},     $totects / $timeweight;
        push @{$outLC{ergrateerr}},  sqrt( $toteerr ) / $timeweight;
        push @{$outLC{bkgrate}},     $totbcts / $timeweight;
        push @{$outLC{bkgrateerr}},  sqrt( $totberr ) / $timeweight;
        push @{$outLC{bkgbackcorr}}, $totbcor / $timeweight;
        push @{$outLC{fracexp}},     $timeweight / ( $stop - $start );
        push @{$outLC{mincts}},      $cts;
        push @{$outLC{uplimit}},     $upper;
        debug( sprintf( "BIN: %10.6f - %10.6f, %4.2f/%d, %4.2f, %4.2f, %4.2f, %4.2f/%4.2f = %4.2f, %4.2f\n",
                        $start, $stop, $totcts, $cts, sqrt( $toterr ), $totbcts, sqrt( $totberr ), $totbcor, 
                        $timeweight, $totbcor / $timeweight, $snr ) );
    };
    
    #
    # add an upper limit
    #
    my $addupper = sub {
        # set the stop time
        if ( $i - 1 < @{$columns{time}{val}} ) {
            $stop = $columns{time}{val}[ $i - 1 ] + $columns{timedel}{val}[ $i - 1 ] / 2.0;
        } else {
            $stop = $columns{time}{val}[ $i - 2 ] + $columns{timedel}{val}[ $i - 2 ] / 2.0;
        }
   
        if ( $totccor <= 0.0 ) {
            $totccor = 1.0;
        } else {
            $totccor /= $totccwt;
        }
        if ( $totecon <= 0.0 ) {
            $totecon = 5.0e-11;
        } else {
            $totecon /= $totccwt;
        }

        # push the time, rate, etc.
        push @{$outLC{timedel}},     $stop - $start;
        push @{$outLC{time}},        ( $start + $stop ) / 2.0;
        push @{$outLC{rate}},        $upper;
        push @{$outLC{rateerr}},     FLT_NULL;
        push @{$outLC{netrate}},     $upper;
        push @{$outLC{netrateerr}},  FLT_NULL;
        push @{$outLC{corrate}},     $upper * $totccor;
        push @{$outLC{corrateerr}},  FLT_NULL;
        push @{$outLC{ergrate}},     $upper * $totecon;
        push @{$outLC{ergrateerr}},  FLT_NULL;
        push @{$outLC{bkgrate}},     $totbcts / $timeweight;
        push @{$outLC{bkgrateerr}},  sqrt( $totberr ) / $timeweight;
        push @{$outLC{bkgbackcorr}}, $totbcor / $timeweight;
        push @{$outLC{fracexp}},     $timeweight / ( $stop - $start );
        push @{$outLC{mincts}},      $cts;
        push @{$outLC{uplimit}},     $upper;
        debug( sprintf( "UPP: %10.6f - %10.6f, %4.2f/%d, %4.2f, %4.2f, %4.2f, %4.2f/%4.2f = %4.2f, %4.2f\n",
                        $start, $stop, $totcts, $cts, $bgvari, $totbcts, sqrt( $totberr ), $totbcor, 
                        $timeweight, $totbcor / $timeweight, $snr ) );
    };
       
    #
    # merge current counts into last bin - if we are step binning,
    # we assume that merging bins does not drop the S/N below the
    # threshold
    #
    my $mergebin = sub {

        # get the last bin start time, and array pos
        my $j    = $#{$outLC{time}};
        my $strt = $outLC{time}[ $j ] - $outLC{timedel}[ $j ] / 2.0;

        # set the stop time
        if ( $i - 1 < @{$columns{time}{val}} ) {
            $stop = $columns{time}{val}[ $i - 1 ] + $columns{timedel}{val}[ $i - 1 ] / 2.0;
        } else {
            $stop = $columns{time}{val}[ $i - 2 ] + $columns{timedel}{val}[ $i - 2 ] / 2.0;
        }
    
        # add the counts from current bin to last bin
        $totcts     += $oldcts;
        $toterr     += $olderr;
        $totncts    += $oldncts;
        $totnerr    += $oldnerr;
        $totccts    += $oldccts;
        $totcerr    += $oldcerr;
        $totects    += $oldects;
        $toteerr    += $oldeerr;
        $totbcts    += $oldbcts;
        $totberr    += $oldberr;
        $totbcor    += $oldbcor;
        $timeweight += $oldweight;
        
        # replace the last bin
        $outLC{timedel}[ $j ]     = $stop - $strt; 
        $outLC{'time'}[ $j ]      = ( $strt + $stop ) / 2.0;
        $outLC{rate}[ $j ]        = $totcts / $timeweight;
        $outLC{rateerr}[ $j ]     = sqrt( $toterr ) / $timeweight;
        $outLC{netrate}[ $j ]     = $totncts / $timeweight;
        $outLC{netrateerr}[ $j ]  = sqrt( $totnerr ) / $timeweight;
        $outLC{corrate}[ $j ]     = $totccts / $timeweight;
        $outLC{corrateerr}[ $j ]  = sqrt( $totcerr ) / $timeweight;
        $outLC{ergrate}[ $j ]     = $totects / $timeweight;
        $outLC{ergrateerr}[ $j ]  = sqrt( $toteerr ) / $timeweight;
        $outLC{bkgrate}[ $j ]     = $totbcts / $timeweight;
        $outLC{bkgrateerr}[ $j ]  = sqrt( $totberr ) / $timeweight;
        $outLC{bkgbackcorr}[ $j ] = $totbcor / $timeweight;
        $outLC{fracexp}[ $j ]     = $timeweight / ( $stop - $strt );
        $outLC{mincts}[ $j ]      = $cts;
        $outLC{uplimit}[ $j ]     = FLT_NULL;
        debug( sprintf( "MRG: %10.6f - %10.6f, %4.2f/%d, %4.2f, %4.2f, %4.2f, %4.2f/%4.2f = %4.2f, %4.2f\n",
                        $strt, $stop, $totcts, $cts, sqrt( $toterr ), $totbcts, sqrt( $totberr ), $totbcor, 
                        $timeweight, $totbcor / $timeweight, $snr ) );
    };

    #
    # set new start time
    #
    my $setstart = sub {
        if ( $i < @{$columns{time}{val}} ) {
            $start = $columns{time}{val}[ $i ] - $columns{timedel}{val}[ $i ] / 2.0;
        }
    };

    #
    # clear count total and time weight total
    #
    my $cleartotals = sub {
        $totcts  = $toterr = 
        $totncts = $totnerr =
        $totccts = $totcerr =
        $totects = $toteerr =
        $totbcts = $totberr =
        $totccor = $totecon = $totccwt =
        $totbcor = $timeweight = $bgvari = $snr = 0.0;
    };

    #
    # clear saved count total and time weight total
    #
    my $clearold = sub {
        $oldcts  = $olderr = 
        $oldncts = $oldnerr =
        $oldccts = $oldcerr =
        $oldects = $oldeerr =
        $oldbcts = $oldberr =
        $oldbcor = $oldweight = 0.0;
    };

    #
    # save count total
    #
    my $savetotals = sub {
        $oldcts  = $totcts;
        $olderr  = $toterr;
        $oldncts = $totncts;
        $oldnerr = $totnerr;
        $oldccts = $totccts;
        $oldcerr = $totcerr;
        $oldects = $totects;
        $oldeerr = $toteerr;
        $oldbcts = $totbcts;
        $oldberr = $totberr;
        $oldbcor = $totbcor;
        $oldweight = $timeweight;
    };
    # END anonymous subs

    # BIN!
    while ( $row <= $nrows + 1 ) {
       
        if ( $row <= $nrows ) {

            # read the next chunk if we need it
            if ( !$recursing && $i >= @{$columns{time}{val}} ) {
                if ( $row + $optrows > $nrows ) {
                    $optrows = $nrows - $row;
                }
                if ( $optrows > 0 ) {
                    debug( "reading $optrows rows ($row-" . ( $row + $optrows )
                         . ") from file $lcfile\[1\]\n" );
                    foreach my $col ( @colnames ) {
                        # overlap by one row
                        @{$columns{$col}{val}} = ( 0.0 ) x $optrows + 1;
                        $fptr->read_col( $columns{$col}{type}, $columns{$col}{colnum}, $row - 1, 1,
                                         $optrows + 1, $null, \@{$columns{$col}{val}}, $any, $stat );
                    }
                    $i = 1;
                }
            }
            # check for gap
            if ( defined $columns{time}{val}[ $i ] && defined $columns{timedel}{val}[ $i ] ) {
                my $gapstop = $columns{time}{val}[ $i ] - $columns{timedel}{val}[ $i ] / 2.0;
                my $gapstrt = $columns{time}{val}[ $i - 1 ] + $columns{timedel}{val}[ $i - 1 ] / 2.0;
                $hasGap = $gapstop - $gapstrt > $gapintv;
            } else {
                $hasGap = 1;
            }
        } else {
            # always a gap at the end
            $hasGap = 1;
        }

        # if we have a gap, then act accordingly
        if ( $hasGap ) {
          
            # good bin
            if ( $timeweight >= $mintime && 
                      ( ( $step && $totcts >= $cts && $snr >= $detthr ) ||
                        ( !$step && $totcts >= $cts ) ) ) {
                $upper = FLT_NULL;
                $addbin->( );

            # glom onto last bin
            #} elsif ( ( $hasBin && !$recursing && $totcts > 0.0 ) ||
            #          ( @{$outLC{time}} && $recursing && $snr >= $detthr ) ) {
            } elsif ( $hasBin && !$recursing && $totcts > 0.0 ) {
                $mergebin->( );

            # here we are crossing gaps, but only in PC mode
            # and only if recursing
            } elsif ( $recursing && $mode eq 'pc' && $timeweight >= $mintime ) {
                #$clearold->( );
                #$hasBin = 0;
                goto ACCUM;

            # mark as low bin if there are counts
            # or we are not recursing
            } elsif ( ( !$recursing || $totcts > 0.0 ) && $timeweight >= $mintime ) {
                $upper = -1;
                $addbin->( );
            }

            # clear!
            if ( !$recursing ) { $clearold->( ); } else { $savetotals->( ); }
            $cleartotals->( );
            $setstart->( );
            $hasBin = 0;

        # otherwise, just keep binning
        } elsif ( $timeweight >= $mintime && 
                  ( ( $step && $totcts >= $cts && $snr >= $detthr ) ||
                    ( !$step && $totcts >= $cts ) ) ) {
            $upper = FLT_NULL;
            $addbin->( );
            $savetotals->( );
            $cleartotals->( );
            $setstart->( );
            $hasBin = 1;
        }
     
ACCUM:
        # accumulate counts and other quantities
        if ( $i < @{$columns{time}{val}} ) {

            my $scaledTime = $columns{timedel}{val}[ $i ] * $columns{fracexp}{val}[ $i ];
            $timeweight += $scaledTime;
            $totcts  += $columns{rateraw}{val}[ $i ] * $scaledTime;
            $toterr  += ( $columns{errorraw}{val}[ $i ] * $scaledTime )**2.0;
            $totncts += $columns{ratenet}{val}[ $i ] * $scaledTime;
            $totnerr += ( $columns{errornet}{val}[ $i ] * $scaledTime )**2.0;
            $totccts += $columns{ratecor}{val}[ $i ] * $scaledTime;
            $totcerr += ( $columns{errorcor}{val}[ $i ] * $scaledTime )**2.0;
            $totects += $columns{rateflux}{val}[ $i ] * $scaledTime;
            $toteerr += ( $columns{errorflux}{val}[ $i ] * $scaledTime )**2.0;
            $totbcts += $columns{ratebkg}{val}[ $i ] * $scaledTime;
            $totberr += ( $columns{errorbkg}{val}[ $i ] * $scaledTime )**2.0;
            $totbcor += $columns{bkgbackcorr}{val}[ $i ] * $scaledTime;
            if ( $columns{rateraw}{val}[ $i ]  >= 0.0 &&
                 $columns{ratenet}{val}[ $i ]  >  0.0 &&
                 $columns{ratecor}{val}[ $i ]  >= 0.0 &&
                 $columns{rateflux}{val}[ $i ] >= 0.0 ) {
                $totccor += $columns{rateraw}{val}[ $i ] *
                            $columns{ratecor}{val}[ $i ] / $columns{ratenet}{val}[ $i ];
                $totecon += $columns{rateraw}{val}[ $i ] *
                            $columns{rateflux}{val}[ $i ] / $columns{ratenet}{val}[ $i ];
                $totccwt += $columns{rateraw}{val}[ $i ];
            }

            # accum the bkg variance
            my $scrate = $columns{ratebkg}{val}[ $i ] * $scaledTime;
            $bgvari   += $scrate * $columns{bkgbackcorr}{val}[ $i ];

            # get S/N
            $snr = ( $totcts - $totbcts );
            if ( $snr + $bgvari > 0.0 ) {
                $snr /= sqrt( $totcts + $bgvari );
            } else {
                $snr = 0.0;
            }

            # if step binning adjust bin size based on raw rate
            if ( $step && $timeweight > 0 ) {
                my $j = 0;
                my $rate = $totcts / $timeweight;
                while ( $j < @$rates && $rate > $rates->[ $j ] ) { $j++; }
                $cts = $counts->[ $j ];
            }
        }
        $row++;
        if ( $recursing ) {
            debug( "\$i = $i\t\$row = $row\t\$nrows = $nrows\t\$totcts = $totcts\n" );
        }
        $i++;
    }
    debug( "\$i = $i\t\$row = $row\t\$nrows = $nrows\n" );
    debug( "scalar(\@{\$columns{time}{val}}) = " . scalar( @{$columns{time}{val}} ) );

    debug( "checking for left-over data after binning\n" );
    debug( "\$i = $i\t\$row = $row\t\$nrows = $nrows\n" );
    debug( "scalar(\@{\$columns{time}{val}}) = " . scalar( @{$columns{time}{val}} ) );
    
    # set the stop time
    $stop = $columns{time}{val}[ $i - 2 ] + $columns{timedel}{val}[ $i - 2 ] / 2.0;

    debug( "start of left-over data: $start\n" );
    debug( "stop  of left-over data: $stop\n" );
    debug( "total left-over counts:  $totcts\n" );
    debug( "total left-over time:    $timeweight\n" );

    # take care of any remaining time
    if ( $stop - $start > 0.0 && $timeweight > 0.0 ) {
       
        # get S/N
        $snr = ( $totcts - $totbcts );
        if ( $snr + $bgvari > 0.0 ) {
            $snr /= sqrt( $totcts + $bgvari );
        } else {
            $snr = 0.0;
        }

        # adjust bin size based on raw rate
        if ( $step ) {
            my $j = 0;
            my $rate = $totcts / $timeweight;
            while ( $j < @$rates && $rate > $rates->[ $j ] ) { $j++; }
            $cts = $counts->[ $j ];
        }

        # new detection
        if ( $timeweight >= $mintime && 
             ( ( $step && $totcts >= $cts && $snr >= $detthr ) ||
               ( !$step && $totcts >= $cts ) ) ) {
            $upper = FLT_NULL;
            $addbin->( );
        # merged bin
        } elsif ( @{$outLC{time}} && $mode eq 'pc' && $snr >= $detthr ) {
            $mergebin->( );
        # low bin
        } elsif ( $snr >= $detthr ) {
            $upper = -1;
            $addbin->( );
        # upper limit in PC mode or either mode if cutlowbins=no
        } elsif ( $mode eq 'pc' || ( $mode eq 'wt' &&
                                     !$params{cutlowbins} &&
                                     !$params{cutlowbinswt} ) ) {
            # get a new upper limit
            my $stat = 0;
            ( $stat, $upper ) = getUpperLimit( undef, $totcts, $totbcts, 
                                               $detthr );
            return $stat unless $stat == 0;
            $upper /= $timeweight;
            $addupper->( );
        # WT mode && cutlowbins - mark as bad for the cutting
        } else {
            $upper = -1;
            $addbin->( );
        }
    }

    # something went wrong
    return ( $stat, \%outLC ) unless $stat == 0;

    # don't recurse again
    if ( $recursing ) {
        return ( $stat, \%outLC );
    }
    
    # recurse to merge low bins if possible
    # but first get the initial counts/bin if step binnning
    if ( $step ) {
        my $j = 0;
        my $rate = $outLC{rate}[ 0 ];
        while ( $j < @$rates && $rate > $rates->[ $j ] ) { $j++; }
        $cts = $counts->[ $j ];
    }
    my $newLC;
    
    ( $stat, $newLC ) = binLCCts( \%outLC, $cts, $gapintv, $mode,
                                  $mintime, $step, $detthr, $rates, $counts );
    return ( $stat, $newLC );
}

#
# cutLowBins -
#
#       cuts low bins (marked as ->{uplimit} == -1)
#
sub cutLowBins {

    my $raw = shift;

    my $k = scalar( @{$raw->{time}} ) - 1;

    chat( 2, "handling low count bins\n" );
    while ( $k >= 0 ) {
        if ( $raw->{uplimit}[ $k ] == -1  ) {
            chat( 3, "purging low bin at $raw->{time}[ $k ] s after trigger\n" );
            splice @{$raw->{rate}},    $k, 1;
            splice @{$raw->{rateerr}}, $k, 1;
            splice @{$raw->{time}},    $k, 1;
            splice @{$raw->{timedel}}, $k, 1;
            splice @{$raw->{fracexp}}, $k, 1;
            splice @{$raw->{uplimit}}, $k, 1;
            
            splice @{$raw->{bkgrate}},    $k, 1;
            splice @{$raw->{bkgrateerr}}, $k, 1;
            splice @{$raw->{bkgbackcorr}},$k, 1;

            splice @{$raw->{netrate}},    $k, 1;
            splice @{$raw->{netrateerr}}, $k, 1;
            
            splice @{$raw->{corrate}},    $k, 1;
            splice @{$raw->{corrateerr}}, $k, 1;
            
            splice @{$raw->{ergrate}},    $k, 1;
            splice @{$raw->{ergrateerr}}, $k, 1;

            if ( @{$raw->{mincts}} ) {
                splice @{$raw->{mincts}}, $k, 1;
            }
        }
        $k--;
    }
}

#
# applyHardToSoft -
#
#       applies hard curve bin boundaries to soft curve
#
sub applyHardToSoft {

    my ( $hardlc, $softlcfile, $mode ) = @_;
    
    chat( 2, "applying hard band bins to soft band light curve\n" );
 
    # read the start of the soft lc
    my $any = 0;
    my $null = FLT_NULL;
    my $stat = 0;
    my $nrows;
    my $optrows;
    my $fptr = Astro::FITS::CFITSIO::open_file( "$softlcfile\[1\]", READONLY, $stat );
    $fptr->get_num_rows( $nrows, $stat );
    if ( $nrows < 1 ) { 
        warnhi( 1, "no rows in file $softlcfile\[1\]\n" );
        return 0;
    }
    $fptr->get_rowsize( $optrows, $stat );
    if ( $optrows > $nrows ) { $optrows = $nrows; }

    my %columns = ( );
    my @colnames = qw( time timedel fracexp ratenet errornet ratecor 
                       errorcor rateflux errorflux rateraw errorraw ratebkg errorbkg bkgbackcorr );
    foreach my $col ( @colnames ) {
        my $num;
        $fptr->get_colnum( CASEINSEN, $col, $num, $stat );
        $columns{$col}{colnum} = $num;
        
        my $type;
        $fptr->get_coltype( $num, $type, undef, undef, $stat );
        $columns{$col}{type} = $type;

        @{$columns{$col}{val}} = ( );
        $fptr->read_col( $columns{$col}{type}, $columns{$col}{colnum}, 1, 1,
                         $optrows, $null, \@{$columns{$col}{val}}, $any, $stat );
    }

    my %softlc = getEmptyLC( );
    
    # loop over hard lc time
    my $i = 0; # hard lc
    my $j = 0; # index of current soft chunk
    my $row = 1; # row number in soft file
    while ( $i < @{$hardlc->{time}} ) {
        my $binstart = $hardlc->{time}[ $i ] - $hardlc->{timedel}[ $i ] / 2.0;
        my $binstop  = $hardlc->{time}[ $i ] + $hardlc->{timedel}[ $i ] / 2.0;
       
        $softlc{time}[ $i ]        = $hardlc->{time}[ $i ];
        $softlc{timedel}[ $i ]     = $hardlc->{timedel}[ $i ];
        $softlc{fracexp}[ $i ]     = $hardlc->{fracexp}[ $i ];
        $softlc{rate}[ $i ]        = 0.0;
        $softlc{rateerr}[ $i ]     = 0.0;
        $softlc{netrate}[ $i ]     = 0.0;
        $softlc{netrateerr}[ $i ]  = 0.0;
        $softlc{corrate}[ $i ]     = 0.0;
        $softlc{corrateerr}[ $i ]  = 0.0;
        $softlc{ergrate}[ $i ]     = 0.0;
        $softlc{ergrateerr}[ $i ]  = 0.0;
        $softlc{bkgrate}[ $i ]     = 0.0;
        $softlc{bkgrateerr}[ $i ]  = 0.0;
        $softlc{bkgbackcorr}[ $i ] = 0.0;
        my $realtime = 0.0;
        my $bgvari   = 0.0;
        
        # accumulate rate over bin span
        my $looped = 0;
        while ( $row <= $nrows ) {
            # read the next chunk if we need it
            if ( $j >= $optrows - 1 ) {
                if ( $row + $optrows > $nrows ) {
                    $optrows = $nrows - $row;
                }
                if ( $optrows != 0 ) {
                    debug( "reading rows $row-" . ( $row + $optrows - 1 )
                         . " from file $softlcfile\[1\]\n" );
                    foreach my $col ( @colnames ) {
                        @{$columns{$col}{val}} = ( 0.0 ) x $optrows;
                        $fptr->read_col( $columns{$col}{type}, $columns{$col}{colnum}, $row, 1,
                                         $optrows, $null, \@{$columns{$col}{val}}, $any, $stat );
                    }
                    $j = 0;
                } else {
                    last;
                }
            }
            if ( $columns{time}{val}[ $j ] > $binstop ) {
                $j++;
                $row++;
                last;
            }
            
            if ( $columns{time}{val}[ $j ] < $binstart ) {
                $j++;
                $row++;
                if ( $row > $nrows ) { last; }
                next; 
            }
            my $scale  = $columns{timedel}{val}[ $j ] * $columns{fracexp}{val}[ $j ];
            $softlc{rate}[ $i ]        += $columns{rateraw}{val}[ $j ] * $scale;
            $softlc{rateerr}[ $i ]     += ( $columns{errorraw}{val}[ $j ] * $scale )**2.0;
            $softlc{netrate}[ $i ]     += $columns{ratenet}{val}[ $j ] * $scale;
            $softlc{netrateerr}[ $i ]  += ( $columns{errornet}{val}[ $j ] * $scale )**2.0;
            $softlc{corrate}[ $i ]     += $columns{ratecor}{val}[ $j ] * $scale;
            $softlc{corrateerr}[ $i ]  += ( $columns{errorcor}{val}[ $j ] * $scale )**2.0;
            $softlc{ergrate}[ $i ]     += $columns{rateflux}{val}[ $j ] * $scale;
            $softlc{ergrateerr}[ $i ]  += ( $columns{errorflux}{val}[ $j ] * $scale )**2.0;
            $softlc{bkgrate}[ $i ]     += $columns{ratebkg}{val}[ $j ] * $scale;
            $softlc{bkgrateerr}[ $i ]  += ( $columns{errorbkg}{val}[ $j ] * $scale )**2.0;
            $softlc{bkgbackcorr}[ $i ] += $columns{bkgbackcorr}{val}[ $j ] * $scale;
            
            # accum the bkg variance
            my $scrate = $columns{ratebkg}{val}[ $j ] * $scale;
            $bgvari   += $scrate * $columns{bkgbackcorr}{val}[ $j ];

            $realtime += $scale;
            $row++;
            $j++;
            $looped++;
        }

        # if we have a bin (which we should) push it
        if ( $looped ) {
            $bgvari /= $realtime;
            my $snr = ( $softlc{rate}[ $i ] - $softlc{bkgrate}[ $i ] );
            if ( $snr + $bgvari > 0.0 ) {
                $snr /= sqrt( $softlc{rate}[ $i ] + $bgvari );
            } else {
                $snr = 0.0;
            }

            # here we consider anything that was detected in the 
            # hard band, to be a detection in the soft band also
            # the hard band should only have an upper limit at the
            # end of the data, so that is the only place the 
            # soft band will either.
            if ( $i < @{$hardlc->{time}} - 1 || $snr >= $params{detthresh} ) {

                # mark as low if hard low
                $softlc{uplimit}[ $i ] = $hardlc->{uplimit}[ $i ];

                # finalize the bin
                debug( sprintf( "BIN: %.14g s - %.14g s, $softlc{rate}[ $i ] cts\n",
                                $hardlc->{time}[ $i ] - $hardlc->{timedel}[ $i ] / 2.0,
                                $hardlc->{time}[ $i ] + $hardlc->{timedel}[ $i ] / 2.0 ) );
                $softlc{rate}[ $i ]        /= $realtime;
                $softlc{rateerr}[ $i ]      = sqrt( $softlc{rateerr}[ $i ] ) / $realtime;
                $softlc{netrate}[ $i ]     /= $realtime;
                $softlc{netrateerr}[ $i ]   = sqrt( $softlc{netrateerr}[ $i ] ) / $realtime;
                $softlc{corrate}[ $i ]     /= $realtime;
                $softlc{corrateerr}[ $i ]   = sqrt( $softlc{corrateerr}[ $i ] ) / $realtime;
                $softlc{ergrate}[ $i ]     /= $realtime;
                $softlc{ergrateerr}[ $i ]   = sqrt( $softlc{ergrateerr}[ $i ] ) / $realtime;
                $softlc{bkgrate}[ $i ]     /= $realtime;
                $softlc{bkgrateerr}[ $i ]   = sqrt( $softlc{bkgrateerr}[ $i ] ) / $realtime;
                $softlc{bkgbackcorr}[ $i ] /= $realtime;
            } else {

                # get upper limit
                my $upper;
                ( $stat, $upper ) = getUpperLimit( undef, $softlc{rate}[ $i ], $softlc{bkgrate}[ $i ], 
                                                   $params{detthresh} );
                return $stat unless $stat == 0;
                $upper /= $realtime;
                debug( sprintf( "UPP: %.14g s - %.14g s, $upper c/s\n",
                                $hardlc->{time}[ $i ] - $hardlc->{timedel}[ $i ] / 2.0,
                                $hardlc->{time}[ $i ] + $hardlc->{timedel}[ $i ] / 2.0 ) );
                $softlc{rate}[ $i ]         = $upper;
                $softlc{rateerr}[ $i ]      = FLT_NULL;
                if ( $softlc{netrate}[ $i ] != 0.0 ) {
                    $softlc{corrate}[ $i ] = $upper * $softlc{corrate}[ $i ] / $softlc{netrate}[ $i ];
                } else {
                    $softlc{corrate}[ $i ] = $upper;
                }
                $softlc{corrateerr}[ $i ]   = FLT_NULL;
                if ( $softlc{netrate}[ $i ] != 0.0 ) {
                    $softlc{ergrate}[ $i ] = $upper * $softlc{ergrate}[ $i ] / $softlc{netrate}[ $i ];
                } else {
                    my $j = $i;
                    # determine the flux conversion from earlier bins
                    while ( $j >= 0 ) {
                        if ( $softlc{netrate}[ $j ] != 0.0 ) {
                            $softlc{ergrate}[ $i ] = $upper * $softlc{ergrate}[ $j ] / $softlc{netrate}[ $j ];
                            last;
                        }
                        $j--;
                    }
                    if ( $j == -1 ) {
                        warnhi( 1, "Could not determine flux conversion for soft band upper limit!\n" );
                        warnhi( 1, "Soft band upper limit is not reliable!\n" );
                        $softlc{ergrate}[ $i ] = $hardlc->{ergrate}[ $i ];
                    }
                }
                $softlc{ergrateerr}[ $i ]   = FLT_NULL;
                $softlc{netrate}[ $i ]      = $upper;
                $softlc{netrateerr}[ $i ]   = FLT_NULL;
                $softlc{bkgrate}[ $i ]     /= $realtime;
                $softlc{bkgrateerr}[ $i ]   = sqrt( $softlc{bkgrateerr}[ $i ] ) / $realtime;
                $softlc{bkgbackcorr}[ $i ] /= $realtime;

                # mark as low if hard low
                $softlc{uplimit}[ $i ] = $hardlc->{uplimit}[ $i ];
            }
        }
        $i++;
    }
    $fptr->close_file( $stat );
    return ( $stat, \%softlc );
}

sub writeXRTGRBLC {
    my ( $lc, $file, $ra, $dec, $trigtime, $lctempl ) = @_;

    my ( $stat, $fptr ) = beginLightCurveFITS( $file, $ra, $dec, $trigtime, $lctempl, 0 );
    return $stat unless $stat == 0;

    $stat = appendLightCurveFITS( $fptr, $lc );
    return $stat unless $stat == 0;

    fits_close_file( $fptr, $stat );
    return $stat
}

sub beginLightCurveFITS {
    my ( $file, $ra, $dec, $trigtime, $lctempl, $addbackcorr ) = @_;
    
    my ( $fptr, $hdutype, $col, $code, $rpt, $wdt, $any );
    
    my $stat = 0;

    my $nrows;
    
    # keywords for primary header and rate table
    my @primkeys = ( 
        { name => 'TELESCOP', type => TSTRING, 
          comm => 'Telescope (mission) name' }, 
        { name => 'INSTRUME', type => TSTRING, 
          comm => 'Intrument name' }, 
        { name => 'OBJECT',   type => TSTRING, val => $params{object},
          comm => 'Name of observed object' },
        { name => 'ORIGIN',   type => TSTRING,
          comm => 'Origin of fits file' }, 
        { name => 'TSTART',   type => TDOUBLE,
          comm => 'Start time' }, 
        { name => 'TSTOP',    type => TDOUBLE,
          comm => 'Stop time' }, 
        { name => 'TELAPSE',  type => TDOUBLE,
          comm => 'Elapsed time' },
        { name => 'MJDREFI',  type => TINT,
          comm => 'MJD reference day' },
        { name => 'MJDREFF',  type => TDOUBLE, 
          comm => 'MJD reference (fraction of day)' },
        { name => 'TIMEREF',  type => TSTRING, 
          comm => 'Reference time' },
        { name => 'TIMESYS',  type => TSTRING, 
          comm => 'Time measured from' },
        { name => 'TIMEUNIT', type => TSTRING, 
          comm => 'Unit of time keywords' }
    );

    # keywords for rate table only (primkeys are also written)
    my @ratekeys = (
        { name => 'TIMEZERO',  type => TDOUBLE,
          comm => 'Time zero' },
        { name => 'TASSIGN',   type => TSTRING,
          comm => 'Time assigned by clock' },
        { name => 'TIERRELA',  type => TDOUBLE,
          comm => '[s/s] relative errors expressed as rate' },
        { name => 'TIERABSO',  type => TDOUBLE,
          comm => '[s] timing precision in seconds' },
        { name => 'DEADAPP',   type => TLOGICAL, val  => 1,
          comm => 'Has DEADC been applied to data' },
        { name => 'HDUCLASS',  type => TSTRING, val  => 'OGIP',
          comm => 'Format conforms to OGIP/GSFC conventions' },
        { name => 'HDUCLAS1',  type => TSTRING, val  => 'LIGHTCURVE',
          comm => 'Extension contains a light curve' },
        { name => 'PROCVER',   type => TSTRING,
          comm => 'Processing script version' },
        { name => 'RA_OBJ',    type => TDOUBLE, val  => $ra,
          comm => '[deg] R.A. Object' },
        { name => 'DEC_OBJ',   type => TDOUBLE, val  => $dec,
          comm => '[deg] Dec Object' },
        { name => 'PHALCUT',   type => TINT, 
          comm => 'Minimum PI channel' },
        { name => 'PHAHCUT',   type => TINT,     
          comm => 'Maximum PI channel' },
        { name => 'DATAMODE',  type => TSTRING,
          comm => 'Datamode' }
    ); 
    if ( $params{usetrigtime} ) {
        if ( $params{trigfrombat} ) {
            push @primkeys, { name => 'TRIGTIME', type => TDOUBLE, val => $trigtime,
                              comm => '[s] MET TRIGger Time for Automatic Target' };
        } else {
            push @primkeys, { name => 'TRIGTIME', type => TDOUBLE, val => $trigtime,
                              comm => '[s] MET TRIGger Time from Other Observatory' };
        }
    }
    if ( $params{usetargid} ) {
        push @primkeys, { name => 'TARG_ID', type => TLONG, comm => 'Target ID' };
    }

    # read the template header
    debug( "reading template file for header keywords: $lctempl\n" );
    $fptr = Astro::FITS::CFITSIO::open_file( $lctempl, READONLY, $stat );
    return $stat unless $stat == 0;

    my $primheader = $fptr->read_header( $stat );
    
    $fptr->movnam_hdu( BINARY_TBL, "RATE", 0, $stat );
    
    my $rateheader = $fptr->read_header( $stat );
    
    $fptr->close_file( $stat );
    
    # setup output columns
    my @ttype = qw( TIME XAX_E TIMEDEL FRACEXP RATENET ERRORNET 
                    RATECOR ERRORCOR RATEFLUX ERRORFLUX RATERAW 
                    ERRORRAW RATEBKG ERRORBKG );
    my @tform = qw( 1D 1D 1D 1E 1E 1E 1E 1E 1E 1E 1E 1E 1E 1E );
    my @tunit =   ( 's', 's', 's', 's', 'count/s', 'count/s', 
                    'count/s', 'count/s', 'erg/cm**2/s', 'erg/cm**2/s',
                    'count/s', 'count/s', 'count/s', 'count/s' );
    if ( $addbackcorr ) {
        push @ttype, "BKGBACKCORR";
        push @tform, "1E",
        push @tunit, '';
    }
    my $nfld  = scalar( @ttype );
    
    # create the lightcurve
    $fptr = undef;
    debug( "creating empty FITS table in $file\n" );
    fits_create_file( $fptr, $file, $stat );
    return $stat unless $stat == 0;
    
    # null primary array
    fits_write_grphdr( $fptr, 1, 8, 0, [ 0 ], 0, 1, 1, $stat );
    
    # write primary header keywords
    debug( "writing keywords to primary header\n" );
    foreach my $key ( @primkeys ) {

        # if we know the val, then write it
        if ( exists $key->{val} ) {
            fits_update_key( $fptr, $key->{type}, $key->{name}, 
                            $key->{val},  $key->{comm}, $stat );
        # otherwise, default to the rate table extension
        } elsif ( exists $rateheader->{$key->{name}} ) {
            $rateheader->{$key->{name}} =~ s/[\"\']//g;
            fits_update_key( $fptr, $key->{type}, $key->{name},
                            $rateheader->{$key->{name}}, $key->{comm}, $stat );
        } elsif ( exists $primheader->{$key->{name}} ) {
            $primheader->{$key->{name}} =~ s/[\"\']//g;
            fits_update_key( $fptr, $key->{type}, $key->{name},
                            $primheader->{$key->{name}}, $key->{comm}, $stat );
        }
    }

    # add the rate table
    debug( "inserting BINTABLE extension\n" );
    fits_insert_btbl( $fptr, 0, $nfld, \@ttype, \@tform, \@tunit, "RATE", 0, $stat );
    
    # add the rate table header keywords
    debug( "writing keywords to BINTABLE extension\n" );
    foreach my $key ( @primkeys, @ratekeys ) {

        # if we know the val, then write it
        if ( exists $key->{val} ) {
            fits_update_key( $fptr, $key->{type}, $key->{name}, 
                            $key->{val},  $key->{comm}, $stat );
        # otherwise, default to the rate table extension
        } elsif ( exists $rateheader->{$key->{name}} ) {
            $rateheader->{$key->{name}} =~ s/[\"\']//g;
            fits_update_key( $fptr, $key->{type}, $key->{name},
                            $rateheader->{$key->{name}}, $key->{comm}, $stat );
        } elsif ( exists $primheader->{$key->{name}} ) {
            $primheader->{$key->{name}} =~ s/[\"\']//g;
            fits_update_key( $fptr, $key->{type}, $key->{name},
                            $primheader->{$key->{name}}, $key->{comm}, $stat );
        }
    }
    return ( $stat, $fptr );
}

sub writeUnbinnedLightCurve {

    my ( $fptr, $srclc, $bkglc, $evt, $intnum, $band ) = @_;

    my ( $fptr2, $fptr3 );
    my $stat = 0;
    my $nrows = 0;
    my $optrows = 0;
    my $any;

    # open source/background light curves
    $fptr2 = Astro::FITS::CFITSIO::open_file( $srclc, READONLY, $stat );
    $fptr2->movnam_hdu( BINARY_TBL, "RATE", 0, $stat );
    $fptr3 = Astro::FITS::CFITSIO::open_file( $bkglc, READONLY, $stat );
    $fptr3->movnam_hdu( BINARY_TBL, "RATE", 0, $stat );

    # here we assume they have the same # of rows (they better!)
    $fptr2->get_num_rows( $nrows, $stat );
    $fptr2->get_rowsize( $optrows, $stat );

    # get the cols 
    my %srccolumns = ( );
    my %bkgcolumns = ( );
    my @colnames = qw( time rate error fracexp );
    foreach my $col ( @colnames ) {
        my $num;
        my $typ;
        
        $fptr2->get_colnum( CASEINSEN, "$col*", $num, $stat );
        $srccolumns{$col}{colnum} = $num;
        
        $fptr2->get_coltype( $num, $typ, undef, undef, $stat );
        $srccolumns{$col}{type} = $typ;
        
        $fptr3->get_colnum( CASEINSEN, "$col*", $num, $stat );
        $bkgcolumns{$col}{colnum} = $num;
        
        $fptr3->get_coltype( $num, $typ, undef, undef, $stat );
        $bkgcolumns{$col}{type} = $typ;
    }

    # here we assume they have the same timedel and tzero (they better!)
    my ( $tdel, $tzero, $deadc, $deadapp );
    $fptr2->read_key_dbl( 'TIMEDEL',  $tdel,  undef, $stat );
    $fptr2->read_key_dbl( 'TIMEZERO', $tzero, undef, $stat );
    $fptr2->read_key_dbl( 'DEADC',    $deadc, undef, $stat );
    $fptr2->read_key_log( 'DEADAPP',  $deadapp, undef, $stat );

    # assuming this comes from DEADAPP, since everything else better be there
    if ( $stat == KEY_NO_EXIST ) {
        fits_clear_errmsg( );
        $deadapp = 0;
        $stat = 0;
    }
    if ( $deadc <= 0.0 ) {
        $deadc = 1.0;
    }

    # invert if > 1 (not likely for swift)
    if ( $deadc > 1.0 ) {
        $deadc = 1.0 / $deadc;
    }

    # get the right energy conversion factor
    my $econv;
    if    ( $band eq 'fb' )   { $econv = $evt->{engconv}[ $intnum ]; }
    elsif ( $band eq 'hard' ) { $econv = $evt->{hengconv}[ $intnum ]; }
    elsif ( $band eq 'soft' ) { $econv = $evt->{sengconv}[ $intnum ]; }

    # calc the rows
    my $row = 1;
    while ( $row <= $nrows ) {
        my %slc = getEmptyLC( );
        my %blc = getEmptyLC( );
        my %rawLC = getEmptyLC( );
        
        # read a chunk
        if ( $row + $optrows > $nrows ) {
            $optrows = $nrows - $row + 1;
        }
        foreach my $col ( @colnames ) {
            my $nam = $col eq 'error' ? 'rateerr' : $col;
            @{$slc{$nam}} = ( 0.0 ) x $optrows;
            @{$blc{$nam}} = ( 0.0 ) x $optrows;
            $fptr2->read_col( $srccolumns{$col}{type}, $srccolumns{$col}{colnum}, $row, 1,
                              $optrows, undef, \@{$slc{$nam}}, $any, $stat );
            $fptr3->read_col( $bkgcolumns{$col}{type}, $bkgcolumns{$col}{colnum}, $row, 1,
                              $optrows, undef, \@{$blc{$nam}}, $any, $stat );
        }
        return $stat unless $stat == 0;

        for ( my $j = 0; $j < @{$slc{time}}; $j++ ) {
          
            # correct for dead time
            if ( !$deadapp ) {
                $slc{rate}[ $j ]    /= $deadc;
                $slc{rateerr}[ $j ] /= $deadc;
                $blc{rate}[ $j ]    /= $deadc;
                $blc{rateerr}[ $j ] /= $deadc;
            }

            # calc bin time
            my $time = $slc{time}[ $j ] + $tzero - $params{trigtime};
           
            # corrected background light curve
            my $bkg  = $blc{rate}[ $j ] * $evt->{backcorr}[ $intnum ];
            my $bger = $blc{rateerr}[ $j ] * $evt->{backcorr}[ $intnum ];
            
            # calc the rate and error and push into light curves
            my $rate = $slc{rate}[ $j ] - $bkg;
            
            my $err  = sqrt( $slc{rateerr}[ $j ]**2.0 + $bger**2.0 );
            
  	    if($rate < 0) { 
		# Skip bins with negative rate to avoid a
		# negative correction factor.

		$rate = 0;
		$err = 0;
		$row++;
		next;
	    }

	    push @{$rawLC{bkgrate}},     $bkg;
            push @{$rawLC{bkgrateerr}},  $bger;
            push @{$rawLC{bkgbackcorr}}, $evt->{backcorr}[ $intnum ];

            push @{$rawLC{time}},  $time;
            push @{$rawLC{fracexp}}, $slc{fracexp}[ $j ];
            push @{$rawLC{timedel}}, $tdel;
           
            # raw light curve
            push @{$rawLC{rate}},    $slc{rate}[ $j ];
            push @{$rawLC{rateerr}}, $slc{rateerr}[ $j ];
                    
            # background subtracted LC
            push @{$rawLC{netrate}},    $rate;
            push @{$rawLC{netrateerr}}, $err;
            
            # bgk sub, exposure/PSF corrected LC
	    my ($stat, $psfcorr) = linearlyInterpolate($slc{time}[$j] + $tzero, $evt->{psfcorrtimes}[$intnum], $evt->{psfcorrfacts}[$intnum]);
	    return $stat unless $stat == 0;
            my $corr    = $rate    * $psfcorr * $evt->{expcorr}[ $intnum ];
            my $correrr = $err     * $psfcorr * $evt->{expcorr}[ $intnum ];
            
            push @{$rawLC{corrate}},    $corr;
            push @{$rawLC{corrateerr}}, $correrr;

            # bgk sub, exposure/PSF corrected, flux converted LC
            push @{$rawLC{ergrate}},    $corr    * $econv;
            push @{$rawLC{ergrateerr}}, $correrr * $econv;
            
            # increment the row
            $row++;
        }

        # dump them to file
        if ( @{$slc{time}} ) {
            $stat = appendLightCurveFITS( $fptr, \%rawLC );
            return $stat unless $stat == 0;
        }
    }
    return $stat;
}

sub appendLightCurveFITS {
    
    my ( $fptr, $lc ) = @_;

    my $currow;
    my $stat = 0;
    fits_get_num_rows( $fptr, $currow, $stat );
    $currow++;

    my $nrows = scalar( @{$lc->{time}} );
    if ( $nrows < 1 ) {
        warnhi( 1, "no rows to write in appendLightCurveFITS( )\n" );
        return 0;
    }
    my $fname;
    $fptr->file_name( $fname, $stat );
    debug( "writing $nrows rows of data to $fname\n" );

    # generate more columns
    my @xax_e = ( );
    foreach my $timdel ( @{$lc->{timedel}} ) {
        push @xax_e, $timdel / 2.0;
    }

    # add the columns
    my $col;
    fits_get_colnum( $fptr, CASEINSEN, 'TIME', $col, $stat );
    fits_write_col_dbl( $fptr, $col, $currow, 1, $nrows, \@{$lc->{'time'}}, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'XAX_E', $col, $stat );
    fits_write_col_flt( $fptr, $col, $currow, 1, $nrows, \@xax_e, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'TIMEDEL', $col, $stat );
    fits_write_col_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{timedel}}, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'FRACEXP', $col, $stat );
    fits_write_col_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{fracexp}}, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'RATENET', $col, $stat );
    fits_write_colnull_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{netrate}}, FLT_NULL, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'ERRORNET', $col, $stat );
    fits_write_colnull_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{netrateerr}}, FLT_NULL, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'RATECOR', $col, $stat );
    fits_write_colnull_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{corrate}}, FLT_NULL, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'ERRORCOR', $col, $stat );
    fits_write_colnull_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{corrateerr}}, FLT_NULL, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'RATEFLUX', $col, $stat );
    fits_write_colnull_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{ergrate}}, FLT_NULL, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'ERRORFLUX', $col, $stat );
    fits_write_colnull_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{ergrateerr}}, FLT_NULL, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'RATERAW', $col, $stat );
    fits_write_colnull_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{rate}}, FLT_NULL, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'ERRORRAW', $col, $stat );
    fits_write_colnull_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{rateerr}}, FLT_NULL, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'RATEBKG', $col, $stat );
    fits_write_colnull_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{bkgrate}}, FLT_NULL, $stat );
    fits_get_colnum( $fptr, CASEINSEN, 'ERRORBKG', $col, $stat );
    fits_write_colnull_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{bkgrateerr}}, FLT_NULL, $stat );

    my $nstat = 0;
    fits_get_colnum( $fptr, CASEINSEN, 'BKGBACKCORR', $col, $nstat );
    if ( $nstat == COL_NOT_FOUND ) {
        $stat = $stat == 0 ? 0 : $stat;
        fits_clear_errmsg( );
    } else {
        fits_write_colnull_flt( $fptr, $col, $currow, 1, $nrows, \@{$lc->{bkgbackcorr}}, FLT_NULL, $stat );
    }
    return $stat;
} 

#
# genRatioFITS -
#
#       pastes together the soft and hard band lc's,
#       and ftcalc's the hardness ratio
#
sub genRatioFITS {

    my ( $soft, $hard, $ratio ) = @_;

    my $stat = 0;
           
    # dump col filters to file for "command line shortness"'s sake
    my @softext = ( "TIME;XAX_E;FRACEXP;TIMEDEL;",
                    "RATENET1(1E)=RATENET;#TUNIT#(column units)=\'count/s\';",
                    "ERRORNET1(1E)=ERRORNET;#TUNIT#(column units)=\'count/s\';",
                    "RATECOR1(1E)=RATECOR;#TUNIT#(column units)=\'count/s\';",
                    "ERRORCOR1(1E)=ERRORCOR;#TUNIT#(column units)=\'count/s\';",
                    "RATEFLUX1(1E)=RATEFLUX;#TUNIT#(column units)=\'erg/cm**2/s\';",
                    "ERRORFLUX1(1E)=ERRORFLUX;#TUNIT#(column units)=\'erg/cm**2/s\'" );
    my $softcol = dumpListToTxt( @softext );

    my @hardext = ( "RATENET2(1E)=RATENET;#TUNIT#(column units)=\'count/s\';",
                    "ERRORNET2(1E)=ERRORNET;#TUNIT#(column units)=\'count/s\';",
                    "RATECOR2(1E)=RATECOR;#TUNIT#(column units)=\'count/s\';",
                    "ERRORCOR2(1E)=ERRORCOR;#TUNIT#(column units)=\'count/s\';",
                    "RATEFLUX2(1E)=RATEFLUX;#TUNIT#(column units)=\'erg/cm**2/s\';",
                    "ERRORFLUX2(1E)=ERRORFLUX;#TUNIT#(column units)=\'erg/cm**2/s\'" );
    my $hardcol = dumpListToTxt( @hardext );

    # copy the primary of the hard to new file and append the RATEBIN ext
    # this is to preseve keywords, but not get the unbinned light curve in there
    my $tmphard = getTmpFile( "lc" );
    $stat = ftcopy( "$hard\[0\]", $tmphard, 0 );
    return $stat unless $stat == 0;

    $stat = runSystem( "ftappend \"$hard\[RATEBIN\]\[col \@$hardcol\]\" \"$tmphard\" chatter=1 history=no" );
    return $stat unless $stat == 0;

    # now paste on the soft cols
    my $ftpaste = {
        infile    => "${soft}\[RATEBIN\][col \@$softcol\]",
        pastefile => $tmphard,
        outfile   => $ratio,
        copyall   => 'no',
        clobber   => 'yes',
        chatter   => $params{chatter},
        history   => 'no'
    };
    $stat = runSystem( genCmd( 'ftpaste', $ftpaste ) );
    quit( $stat ) unless $stat == 0;

    # use extended syntax to calc the hardness ratio
    # this should take care of NULLing any upper limit rows
    my @ratext = ( "*;RATIO(1E)=(isnull(ERRORCOR1).OR.isnull(ERRORCOR2))?#null:RATECOR2/RATECOR1;",
                   "#TUNIT#(columns units)=\'ratio\';",
                   "RATIOERROR(1E)=RATIO*sqrt((ERRORCOR2/RATECOR2)**2.+(ERRORCOR1/RATECOR1)**2.);",
                   "#TUNIT#(columns units)=\'ratio\'" );
    my $ratcol = dumpListToTxt( @ratext );

    # copy unto thyself
    $stat = ftcopy( "$ratio\[RATEBIN\]\[col \@$ratcol\]", $ratio, 1 );
    unlink $tmphard, $softcol, $hardcol, $ratcol;
    return $stat;
}
            
#
# addLCComments -
#
#       adds descriptive column descriptions to curves
#
sub addLCComments {

    my ( $lc, $type, $emin, $emax, $espl ) = @_;

    my $editfile = getTmpFile( "txt" );
    
    open EDIT, ">$editfile" or die "failed to open file $editfile\n";
    
    print EDIT <<EOF;
COMMENT
COMMENT Column Descriptions:
COMMENT
COMMENT TIME      [s]          - Time bin center
COMMENT XAX_E     [s]          - Time bin half-width
COMMENT TIMEDEL   [s]          - Time bin full width
COMMENT FRACEXP                - Fractional exposure of time bin
COMMENT
EOF
    
    if ( $type eq 'src' ) {
        my $enstr = sprintf( "(%.3g-%.3g keV)", $emin, $emax );
        
        print EDIT <<EOF;
COMMENT RATENET   [count/s]    $enstr - Bkg. Subtracted Rate
COMMENT ERRORNET  [count/s]    $enstr - Error in RATENET
COMMENT RATECOR   [count/s]    $enstr - Corrected Rate
COMMENT ERRORCOR  [count/s]    $enstr - Error in RATECOR
COMMENT RATEFLUX  [erg/cm^2/s] $enstr - Flux Rate
COMMENT ERRORFLUX [erg/cm^2/s] $enstr - Error in ERRORFLUX
COMMENT RATERAW   [count/s]    $enstr - Raw source rate
COMMENT ERRORRAW  [count/s]    $enstr - Error on RATERAW
COMMENT RATEBKG   [count/s]    $enstr - Scaled/Corrected Bkg. Rate
COMMENT ERRORBKG  [count/s]    $enstr - Error in RATEBKG
EOF
    
    } elsif ( $type eq 'ratio' ) {
        my $enstr1 = sprintf( "(%.3g-%.3g keV)", $emin, $espl );
        my $enstr2 = sprintf( "(%.3g-%.3g keV)", $espl, $emax );
        my $spc1   = " " x ( length( $enstr2 ) - length( $enstr1 ) );
        my $spc2   = " " x ( length( $enstr1 ) - length( $enstr2 ) );
        my $spc3   = " " x length( "$enstr1$spc1" );
        
        print EDIT <<EOF;
COMMENT RATENET1   [count/s]    $enstr1$spc1 - Bkg. Subtracted Rate
COMMENT ERRORNET1  [count/s]    $enstr1$spc1 - Error in RATENET
COMMENT RATECOR1   [count/s]    $enstr1$spc1 - Corrected Rate
COMMENT ERRORCOR1  [count/s]    $enstr1$spc1 - Error in RATECOR
COMMENT RATEFLUX1  [erg/cm^2/s] $enstr1$spc1 - Flux Rate
COMMENT ERRORFLUX1 [erg/cm^2/s] $enstr1$spc1 - Error in ERRORFLUX
COMMENT RATERAW1   [count/s]    $enstr1$spc1 - Raw source rate
COMMENT ERRORRAW1  [count/s]    $enstr1$spc1 - Error on RATERAW
COMMENT
COMMENT RATENET2   [count/s]    $enstr2$spc2 - Bkg. Subtracted Rate
COMMENT ERRORNET2  [count/s]    $enstr2$spc2 - Error in RATENET
COMMENT RATECOR2   [count/s]    $enstr2$spc2 - Corrected Rate
COMMENT ERRORCOR2  [count/s]    $enstr2$spc2 - Error in RATECOR
COMMENT RATEFLUX2  [erg/cm^2/s] $enstr2$spc2 - Flux Rate
COMMENT ERRORFLUX2 [erg/cm^2/s] $enstr2$spc2 - Error in ERRORFLUX
COMMENT RATERAW2   [count/s]    $enstr2$spc2 - Raw source rate
COMMENT ERRORRAW2  [count/s]    $enstr2$spc2 - Error on RATERAW
COMMENT
COMMENT RATIO      [ ]          $spc3 - RATECOR2/RATECOR1
COMMENT ERRORRATIO [ ]          $spc3 - Error in RATIO
EOF
        
    }

    print EDIT <<EOF;
COMMENT
COMMENT NOTE: Upper limits are indicated by an INDEF in all
COMMENT       of the ERROR* columns. In this case the RATE*
COMMENT       columns represent the $params{detthresh}-sigma upper limit
COMMENT
EOF
    
    close EDIT;

    my %fthedit = (
        infile     => $lc,
        keyword    => "\@$editfile",
        operation  => "add",
        value      => 0,
        protect    => 'yes',
        longstring => 'no',
        chatter    => 0,
        history    => 'no'
    );
    my $cmd  = genCmd( 'fthedit', \%fthedit );
    my $stat = runSystem( $cmd );
    unlink $editfile;

    return $stat;
}

#
# updateHeaders -
#
#       updates various header keywords
#
sub updateHeaders {

    my $fix = shift;

    my $stat = 0;

    # get the swift software version
    #my $swver = `swiftversion`;
    #chomp $swver;

    # keywords for primary headers and all tables
    my @primkeys = (
        { name => 'TELESCOP', type => TSTRING, 
          comm => 'Telescope (mission) name' }, 
        { name => 'INSTRUME', type => TSTRING, 
          comm => 'Intrument name' }, 
        { name => 'OBJECT',   type => TSTRING, val => $params{object},
          comm => 'Name of observed object' },
        { name => 'ORIGIN',   type => TSTRING,
          comm => 'Origin of fits file' }, 
        { name => 'PROCVER',   type => TSTRING,
          comm => 'Processing script version' },
        { name => 'RA_OBJ',    type => TDOUBLE, val  => $params{srcra},
          comm => '[deg] R.A. Object' },
        { name => 'DEC_OBJ',   type => TDOUBLE, val  => $params{srcdec},
          comm => '[deg] Dec Object' },
        { name => 'MJDREFI',  type => TINT,
          comm => 'MJD reference day' },
        { name => 'MJDREFF',  type => TDOUBLE, 
          comm => 'MJD reference (fraction of day)' },
        { name => 'TIMEREF',  type => TSTRING, 
          comm => 'Reference time' },
        { name => 'TIMESYS',  type => TSTRING, 
          comm => 'Time measured from' },
        { name => 'TIMEUNIT', type => TSTRING, 
          comm => 'Unit of time keywords' },
        { name => 'DATE-OBS', type => TSTRING,
          comm => 'Start date of observations' },
        { name => 'DATE-END', type => TSTRING,
          comm => 'End date of observations' },
        #{ name => 'SWIFTVER', type => TSTRING, val => $swver,
        #  comm => 'Swift Software Version' }
    );
    if ( $params{usetrigtime} ) {
        if ( $params{trigfrombat} ) {
            push @primkeys, { name => 'TRIGTIME', type => TDOUBLE, val => $params{trigtime},
                              comm => '[s] MET TRIGger Time for Automatic Target' };
        } else {
            push @primkeys, { name => 'TRIGTIME', type => TDOUBLE, val => $params{trigtime},
                              comm => '[s] MET TRIGger Time from Other Observatory' };
        }
    }
    if ( $params{usetargid} ) {
        push @primkeys, { name => 'TARG_ID', type => TLONG,
                          comm => 'Target ID' };
    }

    # timing keys
    my @timekeys = ( 
        { name => 'TSTART',   type => TDOUBLE,
          comm => 'Start time' }, 
        { name => 'TSTOP',    type => TDOUBLE,
          comm => 'Stop time' }, 
        { name => 'TELAPSE',  type => TDOUBLE,
          comm => 'Elapsed time' }
    );

    # keywords for rate table only (primkeys are also written)
    my @ratekeys = (
        { name => 'TIMEZERO',  type => TDOUBLE, val => $params{trigtime},
          comm => 'Time zero' },
        { name => 'TASSIGN',   type => TSTRING,
          comm => 'Time assigned by clock' },
        { name => 'TIERRELA',  type => TDOUBLE,
          comm => '[s/s] relative errors expressed as rate' },
        { name => 'TIERABSO',  type => TDOUBLE,
          comm => '[s] timing precision in seconds' },
        { name => 'DEADAPP',   type => TLOGICAL, val  => 1,
          comm => 'Has DEADC been applied to data' },
        { name => 'HDUCLASS',  type => TSTRING, val  => 'OGIP',
          comm => 'Format conforms to OGIP/GSFC conventions' },
        { name => 'HDUCLAS1',  type => TSTRING, val  => 'LIGHTCURVE',
          comm => 'Extension contains a light curve' },
        { name => 'PHALCUT',   type => TINT, 
          comm => 'Minimum PI channel' },
        { name => 'PHAHCUT',   type => TINT,     
          comm => 'Maximum PI channel' },
        { name => 'DATAMODE',  type => TSTRING,
          comm => 'Datamode' }
    );

    # keywords for GTI only
    my @gtikeys = (
        { name => 'HDUCLASS',  type => TSTRING, val  => 'OGIP',
          comm => 'Format conforms to OGIP/GSFC conventions' },
        { name => 'HDUCLAS1',  type => TSTRING, val  => 'GTI',
          comm => 'File contains Good Time Intervals' },
        { name => 'HDUCLAS2',  type => TSTRING, val  => 'STANDARD',
          comm => 'File contains Good Time Intervals' }
    );
  
    # get the "file creation date"
    my ( $datestring, $timeref );
    Astro::FITS::CFITSIO::fits_get_system_time( $datestring, $timeref, $stat );

    # figure the overall TSTART/TSTOP/TELAPSE from the GTI
    # and get the UTC time string for DATE-OBS and DATE-END
    my %gtis = ( );
    foreach my $file ( keys %{$fix} ) {

        my $gti = exists $fix->{$file}{gti}   ? $fix->{$file}{gti} : undef;
        
        if ( defined $gti && !exists $gtis{$gti} ) {
            chat( 3, "determining TSTART/TSTOP/TELAPSE from GTI file $gti\n" );
            my $nrows;
            my $fptr3 = Astro::FITS::CFITSIO::open_file( $gti, READONLY, $stat );
            return $stat unless $stat == 0;
            $fptr3->movnam_hdu( BINARY_TBL, "STDGTI", 0, $stat );
            $fptr3->get_num_rows( $nrows, $stat );

            my @ary = ( );
            my ( $any, $col, $type );
            $fptr3->get_colnum( CASEINSEN, "START", $col, $stat );
            $fptr3->get_coltype( $col, $type, undef, undef, $stat );
            $fptr3->read_col( $type, $col, 1, 1, 1, undef, \@ary, $any, $stat );
            $gtis{$gti}{TSTART} = $ary[ 0 ];
            $fptr3->get_colnum( CASEINSEN, "STOP", $col, $stat );
            $fptr3->get_coltype( $col, $type, undef, undef, $stat );
            $fptr3->read_col( $type, $col, $nrows, 1, 1, undef, \@ary, $any, $stat );
            $gtis{$gti}{TSTOP} = $ary[ 0 ];
            $fptr3->close_file( $stat );

            $gtis{$gti}{TELAPSE} = $gtis{$gti}{TSTOP} - $gtis{$gti}{TSTART};

            # now get the DATE-OBS/DATE-END
            $gtis{$gti}{'DATE-OBS'} = swifttime( $gtis{$gti}{TSTART}, 'm', 's', 'u', 't' );
            $gtis{$gti}{'DATE-END'} = swifttime( $gtis{$gti}{TSTOP}, 'm', 's', 'u', 't' );
            if ( $gtis{$gti}{'DATE-OBS'} == FLT_NULL ||
                 $gtis{$gti}{'DATE-END'} == FLT_NULL ) {
                return -1;
            }
        }
    }
    
    # find the GTI with the earliest start and the one with the latest stop
    my %gtivals = ( );
    $gtivals{TSTART} = 1.0e20;
    $gtivals{TSTOP}  = -1.0e20;
    foreach my $gti ( keys %gtis ) {
        if ( exists $gtis{$gti} ) {
            if ( $gtis{$gti}{TSTART} < $gtivals{TSTART} ) {
                $gtivals{TSTART}     = $gtis{$gti}{TSTART};
                $gtivals{'DATE-OBS'} = $gtis{$gti}{'DATE-OBS'};
            }
            if ( $gtis{$gti}{TSTOP} > $gtivals{TSTOP} ) {
                $gtivals{TSTOP}      = $gtis{$gti}{TSTOP};
                $gtivals{'DATE-END'} = $gtis{$gti}{'DATE-END'};
            }
        }
    }
    $gtivals{TELAPSE} = $gtivals{TSTOP} - $gtivals{TSTART};

    # store the headers of each template so we don't
    # read them over and over again
    my %primary = ( );
    my %ratehdr = ( );
    chat( 2, "updating header keywords in output files\n" );
    foreach my $file ( keys %{$fix} ) {
        my ( $fptr, $fptr2, $fptr3, $nhdus, $extname );
        chat( 2, "working on $file\n" );

        my $lctempl = exists $fix->{$file}{templ} ? $fix->{$file}{templ} : undef;

        # if we have not already, read the headers from the template
        if ( defined $lctempl && !exists $primary{$lctempl} ) {
            chat( 3, "reading keys from template light curve $lctempl\n" );
            $fptr2 = Astro::FITS::CFITSIO::open_file( $lctempl, READONLY, $stat );
            return $stat unless $stat == 0;
            $fptr2->movabs_hdu( 1, ANY_HDU, $stat );
            $primary{$lctempl} = $fptr2->read_header( $stat );
            $fptr2->movnam_hdu( BINARY_TBL, "RATE", 0, $stat );
            $ratehdr{$lctempl} = $fptr2->read_header( $stat );
            $fptr2->close_file( $stat );
        }

        # apply the keys
        chat( 3, "applying keywords to $file\n" );
        $fptr = Astro::FITS::CFITSIO::open_file( $file, READWRITE, $stat );
        return $stat unless $stat == 0;
        $fptr->get_num_hdus( $nhdus, $stat );
        foreach my $i ( 1..$nhdus ) {
            $fptr->movabs_hdu( $i, ANY_HDU, $stat );
            $fptr->read_key_str( 'EXTNAME', $extname, undef, $stat );
            if ( $stat == KEY_NO_EXIST ) {
                fits_clear_errmsg( );
                $stat = 0;
            }

            # update the date
            $fptr->update_key( TSTRING, 'DATE', $datestring, 'File creation date', $stat );
            
            # determine the keys to write based on extname
            my @keys = ( );
            if ( $extname =~ /XRTINFO/ || $extname =~ /XRTFLD/ ) {
                @keys = @primkeys;
            } elsif ( $extname =~ /RATE/ ) {
                @keys = ( @primkeys, @timekeys, @ratekeys );
            } elsif ( $extname =~ /GTI/ ) {
                # kill the MJDREF, since we'll be using the MJDREF[IF] convention
                $fptr->delete_key( "MJDREF", $stat );
                @keys = ( @primkeys, @timekeys, @gtikeys );
            } else {
                @keys = ( @primkeys, @timekeys );
            }
            foreach my $key ( @keys ) {

                # if we know the val, then write it
                if ( exists $key->{val} ) {
                    debug( "writing $key->{name} = $key->{val}\n" );
                    $fptr->update_key( $key->{type}, $key->{name}, 
                                       $key->{val},  $key->{comm}, $stat );
                # next check the GTI
                } elsif ( exists $gtivals{$key->{name}} ) {
                    debug( "writing $key->{name} = $gtivals{$key->{name}}\n" );
                    $fptr->update_key( $key->{type}, $key->{name},
                                       $gtivals{$key->{name}},
                                       $key->{comm}, $stat );
                # otherwise, default to the template rate table extension
                } elsif ( exists $ratehdr{$lctempl}->{$key->{name}} ) {
                    $ratehdr{$lctempl}->{$key->{name}} =~ s/[\"\']//g;
                    debug( "writing $key->{name} = $ratehdr{$lctempl}->{$key->{name}}\n" );
                    $fptr->update_key( $key->{type}, $key->{name},
                                       $ratehdr{$lctempl}->{$key->{name}},
                                       $key->{comm}, $stat );
                # and lastly the primary template header
                } elsif ( exists $primary{$lctempl}->{$key->{name}} ) {
                    $primary{$lctempl}->{$key->{name}} =~ s/[\"\']//g;
                    debug( "writing $key->{name} = $primary{$lctempl}->{$key->{name}}\n" );
                    $fptr->update_key( $key->{type}, $key->{name},
                                       $primary{$lctempl}->{$key->{name}},
                                       $key->{comm}, $stat );
                }
            }

            # update the checksum
            $fptr->write_chksum( $stat );

            # stamp the history if in primary header
            if ( $i == 1 ) {
                HDpar_stamp( $fptr, 1, $stat );
            }
        }
        $fptr->close_file( $stat );
    }
    return $stat;
}

#
# ftCreate -
#
#       runs the ftcreate tool
#
sub ftCreate {

    my ( $cdf, $dat, $outf, $extname, $hdr ) = @_;

    my $ftcreate = {
        cdfile      => $cdf,
        datafile    => $dat,
        outfile     => $outf,
        tabtyp      => 'BINARY',
        nskip       => 0,
        nrows       => 0,
        morehdr     => 0,
        extname     => $extname,
        anull       => 0,
        inull       => 0,
        clobber     => 'yes',
        chatter     => $params{chatter},
        history     => 'no'
    };
    $ftcreate->{headfile} = $hdr if $hdr;
    my $cmd  = genCmd( 'ftcreate', $ftcreate );
    my $stat = runSystem( $cmd );
    return $stat;
}

#
# fStatistic -
#
#       runs the fstatistic ftool and parses output
#
sub fStatistic {
    my ( $file, $coln, $rows, $min, $max ) = @_;
    $rows = '-'     unless defined $rows;
    $min  = 'INDEF' unless defined $min;
    $max  = 'INDEF' unless defined $max;
    my $fstatistic = {
        infile  => $file,
        colname => $coln,
        rows    => $rows,
        outfile => 'STDOUT',
        maxval  => $min,
        minval  => $max,
        clobber => 'no'
    };

    my $cmd = genCmd( 'fstatistic', $fstatistic );
    my ( $stat, $out );
    if ( $params{chatter} >= 4 ) {
        ( $stat, $out ) = runSystem( $cmd );
    } else {
        ( $stat, $out ) = runSystemNoChat( $cmd );
    }

    my %patts = (
        sum   => qr/^ The sum of the selected column is\s+(\S+)$/,
        mean  => qr/^ The mean of the selected column is\s+(\S+)$/,
        sigma => qr/^ The standard deviation of the selected column is\s+(\S+)$/,
        min   => qr/^ The minimum of selected column is\s+(\S+)$/,
        max   => qr/^ The maximum of selected column is\s+(\S+)$/,
        numb  => qr/^ The number of points used in calculation is\s+(\S+)$/
    );

    my %stats = ( );
    foreach my $line ( @$out ) {
        foreach my $key ( keys %patts ) {
            if ( $line =~ /$patts{$key}/ ) {
                # carefully, now...
                my @ret = eval { 
                    my @msg = ( );
                    $SIG{__WARN__} = sub { @msg = @_; };
                    $stats{$key} = $1 * 1.0;
                    return @msg;
                };
                if ( @ret ) {
                    warn @ret;
                    warnlo( 1, "there were non-numeric values ",
                               "encountered in fstatistic output!\n" );
                    $stats{$key} = 0.0;
                }
                delete $patts{$key};
            }
        }
    }
    return ( $stat, %stats );
}

# linearlyInterpolate($x0, \@x, \@y)
#
#    Returns a linearly interpolated $y0 value at input value $x0
#    given an array of @x and @y values.  When $x0 is outside the
#    range of @x, the nearest neighbor y value is returned.  Assumes
#    that @x is sorted in increasing values.
sub linearlyInterpolate {
    
    my ($x0, $x, $y) = @_; # Grab the x and y lists as references and grab the x0 value.
    my $stat = 0;  # Set to 1 if a problem occurs later.
    my $y0 = 0; # Initialize the output y value.
    my $i0 = 0; # The @x array index.

    my $n_x = scalar(@{$x}); # Number of array elements
    my $n_y = scalar(@{$y}); # Number of array elements

    if($n_x <= 0 || $n_y <= 0 || $n_x != $n_y) {
	# Flag a failure if the dimensions of the input @$x or @$y arrays are zero
	# or aren't equal.
	$stat = 1;
    } elsif($n_x == 1) {
	# If there is only one point, return the y value.
	$i0 = 0;
	$y0 = $$y[$i0];
    } elsif($x0 <= $$x[0]) {
	# If $x0 is less than the min of @$x, return the zeroth value of @$y.
	$i0 = 0;
	$y0 = $$y[$i0];
    } elsif($x0 >= $$x[$n_x - 1]) {
	# If $x0 is greater than the max of @$x, return the last value of @$y.
	$i0 = $n_x - 1;
	return ($stat, $$y[$i0]);
    } else {
	# Linearly interpolate.
	# Find i0 so that x[i0] <= x0 < x[i0 + 1].
	$i0 = 0;
	while($i0 < $n_x - 1 && $$x[$i0 + 1] < $x0) {
	    $i0++;
	}
	# Calculate the the interpolated y value $y0.
	$y0 = $$y[$i0] + ($x0 - $$x[$i0])/($$x[$i0 + 1] - $$x[$i0]) * ($$y[$i0 + 1] - $$y[$i0]);
    }
    
    #A few debug lines for checking the interpolation.
    #print "  x[$i0] = $$x[$i0]    y[$i0] = $$y[$i0]\n";
    #print "    x0 = $x0           y0 = $y0\n";
    #if($i0 + 1 < $n_x) {print "  x[", $i0 + 1, "] = ", $$x[$i0 + 1], "    y[", $i0 + 1, "] = ", $$y[$i0 + 1], "\n";}

    return ($stat, $y0);
}

# getMean
# 
#       Returns the mean of an array.
#
sub getMean {
    my $array = $_[0];
    my $sum = 0;
    my $length = scalar(@$array);

    if($length == 0) {
	return undef;
    }

    for(my $i = 0; $i < $length; $i++) {
	$sum += $$array[$i];
    }
    
    return $sum/$length;
}

#
# updateChecksum -
#
#       updates the checksums and datasums of input file
#
sub updateChecksum {

    my $lc = shift;

    # update the checksum and datasum
    my %fchecksum = (
        infile => $lc,
        update => 'yes',
        datasum => 'yes'
    );
    my $cmd = genCmd( 'fchecksum', \%fchecksum );
    my $stat = runSystem( $cmd );
    return $stat;
}

#
# genLCurveQDP -
#
#       generates qdp file and pco file for light curves
#       ! Assumes wt mode comes first !
#
sub genLCurveQDP {

    my ( $fn, $mode, $lc, $emin, $emax ) = @_;

    # plot command file
    my $pco = $fn;
    $pco =~ s/\.qdp$/.pco/;
   
    my @labels = ( 'Rate (count s\\U-1\\D)', 'Corrected Rate (count s\\U-1\\D)', 
                   sprintf( "Flux %3g-%.3g keV (erg cm\\U-2\\d s\\U-1\\D)", $emin, $emax ),
                   'FRACEXP' );
    
    # append to qdp file if not extant
    my $startGroup = 1;
    if ( -e $fn ) {
        open FILE, "<$fn" or die "failed to open file $fn\n";
        my $datagrps = 0;
        while ( my $line = <FILE> ) {
            chomp $line;
            $line =~ s/^\s+//;
            $line =~ s/\s+$//;
            next unless length $line;
            if ( $line =~ /^[+\-0-9\.]/ ) {
                $datagrps++;
            }
            if ( $datagrps > 0 && $line =~ /NO\s+NO\s+NO\s+/ ) {
                $datagrps = 0;
                $startGroup += 4;
            }
        }
        close FILE;
        open FILE, ">>$fn" or die "failed to open $fn for writing\n";
        open PCO, ">>$pco" or die "failed to open $pco for writing\n";
        print  FILE "! " . uc( $mode ) . " mode data\n";
    } else {
        open FILE, ">$fn" or die "failed to open $fn for writing\n";
        open PCO, ">$pco" or die "failed to open $pco for writing\n";

        # header
        print  FILE "! Net, Corrected and Flux light curves produced "
                  . "by $taskName V$taskVersion\n";
        print  FILE "read serr 1 2 3 4\n";
        print  FILE "skip single\n";
        print  FILE "\@" . basename( $pco ) . "\n";
        print  FILE "! " . uc( $mode ) . " mode data\n";
        
        printf FILE "! %10s %12s %12s %12s %12s %12s %12s %12s %12s\n",
                    'Time Since', 'Time Error',
                    'Net Rate',   'Net Rate',
                    'Corr. Rate', 'Corr. Rate',
                    'Flux',       'Flux',
                    'Fractional';
                    
        printf FILE "! %10s %12s %12s %12s %12s %12s %12s %12s %12s\n",
                    'Trigger', '', '', 'Error', '', 'Error', '', 'Error', 'Exposure';
                    
        printf FILE "! %10s %12s %12s %12s %12s %12s %12s %12s %12s\n",
                    '[s]', '[s]', '[count/s]', '[count/s]',
                    '[count/s]', '[count/s]',
                    '[erg/cm^2/s]', '[erg/cm^2/s]', '';
    }
        
    my @uplims = ( );
    my $nodata = 1;
    for ( my $i = 0; $i < @{$lc->{time}}; $i++ ) {
        if ( $lc->{netrateerr}[ $i ] != FLT_NULL ) {
            $nodata = 0;
            printf FILE "%12.6g %12.6g %12.6g %12.6g %12.6g %12.6g %12.6g %12.6g %12.6g\n",
                        $lc->{time}[ $i ], $lc->{timedel}[ $i ] / 2.0,
                        $lc->{netrate}[ $i ], $lc->{netrateerr}[ $i ],
                        $lc->{corrate}[ $i ], $lc->{corrateerr}[ $i ],
                        $lc->{ergrate}[ $i ], $lc->{ergrateerr}[ $i ],
                        $lc->{fracexp}[ $i ];
        } else {

            # save upper limits for separate plot group
            my $str = sprintf "%12.6g %12.6g %12.6g %12s %12.6g %12s %12.6g %12s %12.6g\n",
                        $lc->{time}[ $i ], $lc->{timedel}[ $i ] / 2.0,
                        $lc->{netrate}[ $i ], 'NO',
                        $lc->{corrate}[ $i ], 'NO',
                        $lc->{ergrate}[ $i ], 'NO',
                        $lc->{fracexp}[ $i ];
            push @uplims, $str;
        }
    }
    printf FILE "%12s %12s %12s %12s %12s %12s %12s %12s %12s\n", ( 'NO' ) x 9;
    if ( @uplims ) {
        print FILE "! " . uc( $mode ) . " upper limits\n";
        foreach my $str ( @uplims ) {
            print FILE $str;
        }
        printf FILE "%12s %12s %12s %12s %12s %12s %12s %12s %12s\n", ( 'NO' ) x 9;
    }

    # labels and colors
    if ( !$nodata ) {
        my @indices = ( $startGroup, $startGroup + 1, $startGroup + 2, $startGroup + 3 );
        if ( $mode eq 'pc' ) {
            print PCO "color 2 on " . join ' ', @indices, "\n";
        } else {
            print PCO "color 4 on " . join ' ', @indices, "\n";
        }
        foreach my $i ( 0..3 ) {
            print PCO "lab G$indices[ $i ] $labels[ $i ]\n";
        }
        if ( @uplims ) {
            my @upindices = map( $_ + 4, @indices );
            if ( $mode eq 'pc' ) {
                print PCO "color 2 on " . join ' ', @upindices, "\n";
            } else {
                print PCO "color 4 on " . join ' ', @upindices, "\n";
            }
            print PCO "mark 31 on " . join ' ', @upindices[ 0..2 ], "\n";
            print PCO "mark size 2 on " . join ' ', @upindices[ 0..2 ], "\n";
        }
    } else {
        if ( @uplims ) {
            my @upindices = ( $startGroup, $startGroup + 1, $startGroup + 2, $startGroup + 3 );
            if ( $mode eq 'pc' ) {
                print PCO "color 2 on " . join ' ', @upindices, "\n";
            } else {
                print PCO "color 4 on " . join ' ', @upindices, "\n";
            }
            print PCO "mark 31 on " . join ' ', @upindices[ 0..2 ], "\n";
            print PCO "mark size 2 on " . join ' ', @upindices[ 0..2 ], "\n";
            foreach my $i ( 0..3 ) {
                print PCO "lab G$upindices[ $i ] $labels[ $i ]\n";
            }
        }
    }
    close FILE;
    close PCO;
}

#
# plotCurves -
#
#       runs qdp and plots both corr and flux curves
#
sub plotCurves {

    my $qdp = shift;

    my $qdpbase = basename( $qdp );
   
    # change to plot dir
    my $cwd = getcwd;
    chdir dirname( $qdp );

    my %plotgroups = (
        Corrected => {
            ext => "cb_lc.$params{plotfext}",
            grp => '2 6 10 14'
        },
        Flux      => {
            ext => "fb_lc.$params{plotfext}",
            grp => '3 7 11 15'
        }
    );
    
    foreach my $type ( sort keys %plotgroups ) {
        
        my $plotfile = $qdpbase;
        $plotfile =~ s/ab\.qdp$/$plotgroups{$type}{ext}/;
       
        open QDP, "|qdp $qdpbase 2>&1 > /dev/null" 
            or die "failed to run QDP on $qdp\n";
        print QDP "/null\n";
        
        # black on white
        print QDP "scr white\n";

        # turn off everything to start
        print QDP "yplot off 1 2 3 4 5 6 7 8 9 10 11 12\n";
        
        # turn off annoying things
        print QDP "label file\ntime off\n";
        
        # top window soft rate
        print QDP "win 1\nlocation 0 0.1 1 0.9\n";
        print QDP "yplot on $plotgroups{$type}{grp}\n";
        print QDP "log\ngap error\nres\n";
        if ( $params{usetrigtime} ) {
            if ( $params{trigfrombat} ) {
                print QDP "lab X Time Since BAT Trigger (s)\n";
            } else {
                print QDP "lab X Time Since Trigger (s)\n";
            }
        } else {
            print QDP "lab X Time Since MET $params{trigtime} s (s)\n";
        }

        # title
        print QDP "lab T $type Light Curve\n";

        # plot
        print QDP "cpd ${plotfile}$params{plotftype}\n";
        print QDP "plot\n";
        print QDP "quit\n";

        close QDP;
    }
    chdir $cwd;
}

sub genRatioQDP {

    my ( $fn, $mode, $ratiolc, 
         $emin, $emax, $espl ) = @_;

    # plot command file
    my $pco = $fn;
    $pco =~ s/\.qdp$/.pco/;
  
    # indices for plot groups
    my $softstr   = sprintf "%.3g-%.3g keV", $emin, $espl;
    my $hardstr   = sprintf "%.3g-%.3g keV", $espl, $emax;
    my @labels    = ( "$softstr (ct s\\U-1\\D)", "$hardstr (ct s\\U-1\\D)", 
                      'Ratio', 'FRACEXP' );
  
    # open qdp and pco files, and print header
    my $startGroup = 1;
    if ( -e $fn ) {
        open FILE, "<$fn" or die "failed to open file $fn\n";
        my $datagrps = 0;
        while ( my $line = <FILE> ) {
            chomp $line;
            $line =~ s/^\s+//;
            $line =~ s/\s+$//;
            next unless length $line;
            if ( $line =~ /^[+\-0-9\.]/ ) {
                $datagrps++;
            }
            if ( $datagrps > 0 && $line =~ /NO\s+NO\s+NO\s+/ ) {
                $datagrps = 0;
                $startGroup += 4;
            }
        }
        close FILE;
        open FILE, ">>$fn" or die "failed to open $fn for writing\n";
        open PCO, ">>$pco" or die "failed to open $pco for writing\n";
        print  FILE "! " . uc( $mode ) . " mode data\n";
    } else {
        open FILE, ">$fn" or die "failed to open $fn for writing\n";
        open PCO, ">$pco" or die "failed to open $pco for writing\n";
        
        print  FILE "! $softstr, $hardstr Light Curves, and Hardness Ratio\n";
        print  FILE "! produced by $taskName V$taskVersion\n";
        print  FILE "read serr 1 2 3 4\n";
        print  FILE "skip single\n";
        print  FILE "\@" . basename( $pco ) . "\n";
        print  FILE "! " . uc( $mode ) . " mode data\n";
        
        printf FILE "! %11s %13s %13s %13s %13s %13s %13s %13s %13s\n",
                    'Time Since', 'Time Error',
                    'Corr. Rate', 'Corr. Rate',
                    'Corr. Rate', 'Corr. Rate',
                    'Ratio',      'Ratio',
                    'Fractional';
                    
        printf FILE "! %11s %13s %13s %13s %13s %13s %13s %13s %13s\n",
                    'Trigger', '', '', 'Error', '', 'Error', '', 'Error', 'Exposure';
                    
        printf FILE "! %11s %13s %13s %13s %13s %13s %13s %13s %13s\n",
                    '', '', $softstr, $softstr, $hardstr, $hardstr,'', '', '';

        printf FILE "! %11s %13s %13s %13s %13s %13s %13s %13s %13s\n",
                    '[s]', '[s]', '[count/s]', '[count/s]',
                    '[count/s]', '[count/s]', '', '', '';
    }
        
    # ftlist the file with nice formatting
    my @ext = ( "TIME;#TDISP#='G13.6';XAX_E;#TDISP#='G13.6';",
                "RATECOR1;#TDISP#='G13.6';ERRORCOR1;#TDISP#='G13.6';",
                "RATECOR2;#TDISP#='G13.6';ERRORCOR2;#TDISP#='G13.6';",
                "RATIO;#TDISP#='G13.6';RATIOERROR;#TDISP#='G13.6';FRACEXP;#TDISP#='G13.6'" );
    my $colfile = dumpListToTxt( @ext );
    my $file = "$ratiolc\[RATEBIN\]\[col \@$colfile\]";
    my $ftlist = {
        infile    => $file,
        option    => 'T',
        outfile   => '-',
        columns   => '1,2,4,5,6,7,8,9,3',
        rows      => '-',
        vector    => '-',
        separator => ' ',
        rownum    => 'no',
        colheader => 'no'
    };
    open DATA, genCmd( 'ftlist', $ftlist ) . " |" 
        or die "failed to ftlist file $file\n";

    # parse the output to handle upper limits
    my @uplims = ( );
    my $nodata = 1;
    while ( my $line = <DATA> ) {
        chomp $line;
        $line =~ s/^\s*(.*?)\s*$/$1/;
        $line =~ s/(NULL|INDEF)/NO/g;
        my @split = split /\s+/, $line;
        my @upp   = ( );
        my $numnos = 0;
        for ( my $i = 0; $i < @split; $i++ ) {
            if ( $i < 6 && $split[ $i ] =~ /NO/ ) {
                $split[ $i - 1 ] = 'NO';
                $split[ $i ]     = 'NO';
                $upp[ $i - 1 ]   = $split[ $i - 1];
                $upp[ $i ]       = 'NO';
                $numnos++;
            } else {
                $upp[ $i ] = $split[ $i ];
            }
        }
        if ( $numnos < 2 ) {
            $nodata = 0;
            printf FILE "%13s %13s %13s %13s %13s %13s %13s %13s %13s\n", @split;
        } elsif ( $numnos > 0 ) {
            push @uplims, ( sprintf "%13s %13s %13s %13s %13s %13s %13s %13s %13s\n", @upp );
        }

    }
    close DATA;
    unlink $colfile;

    # separator and upper limits
    printf FILE "%13s %13s %13s %13s %13s %13s %13s %13s %13s\n", ( 'NO' ) x 9;
    if ( @uplims ) {
        print FILE "! " . uc( $mode ) . " upper limits\n";
        foreach my $str ( @uplims ) {
            print FILE $str;
        }
        printf FILE "%13s %13s %13s %13s %13s %13s %13s %13s %13s\n", ( 'NO' ) x 9;
    }

    # labels and colors
    if ( !$nodata ) {
        my @indices = ( $startGroup, $startGroup + 1,
                        $startGroup + 2, $startGroup + 3 );
        if ( $mode eq 'pc' ) {
            print PCO "color 2 on " . join ' ', @indices, "\n";
        } else {
            print PCO "color 4 on " . join ' ', @indices, "\n";
        }
        foreach my $i ( 0..3 ) {
            print PCO "lab G$indices[ $i ] $labels[ $i ]\n";
        }
        if ( @uplims ) {
            my @upindices = map( $_ + 4, @indices );
            if ( $mode eq 'pc' ) {
                print PCO "color 2 on " . join ' ', @upindices, "\n";
            } else {
                print PCO "color 4 on " . join ' ', @upindices, "\n";
            }
            print PCO "mark 31 on " . join ' ', @upindices[ 0..1 ], "\n";
            print PCO "mark size 2 on " . join ' ', @upindices[ 0..1 ], "\n";
        }
    } else {
        if ( @uplims ) {
            my @upindices = ( $startGroup, $startGroup + 1, 
                              $startGroup + 2, $startGroup + 3 );
            if ( $mode eq 'pc' ) {
                print PCO "color 2 on " . join ' ', @upindices, "\n";
            } else {
                print PCO "color 4 on " . join ' ', @upindices, "\n";
            }
            print PCO "mark 31 on " . join ' ', @upindices[ 0..1 ], "\n";
            print PCO "mark size 2 on " . join ' ', @upindices[ 0..1 ], "\n";
            foreach my $i ( 0..3 ) {
                print PCO "lab G$upindices[ $i ] $labels[ $i ]\n";
            }
        }
    }
    close FILE;
    close PCO;
}

sub plotRatios {

    my ( $qdp, $min, $max ) = @_;

    my $qdpbase = basename( $qdp );
    
    my $plotfile = $qdpbase;
    $plotfile =~ s/\.qdp$/_lc.$params{plotfext}/;
   
    # change to plot dir
    my $cwd = getcwd;
    chdir dirname( $qdp );

    open QDP, "|qdp $qdpbase 2>&1 > /dev/null" 
        or die "failed to run QDP on $qdp\n";
    print QDP "/null\n";
   
    # black on white
    print QDP "scr white\n";

    # turn off FRACEXP
    print QDP "yplot off 4 8 12\n";
    
    # turn off annoying things
    print QDP "label file\ntime off\n";
    
    # top window soft rate
    print QDP "win 1\nlocation 0 0.6 1 0.933\n";
    print QDP "yplot on 1 5 9\n";
    print QDP "log y\nres y $min $max\nlab NX off\n";

    # title
    print QDP "lab T Hardness Ratio\n";

    # middle window hard rate
    print QDP "win 2\nlocation 0 0.33333 1 0.66667\n";
    print QDP "yplot on 2 6 10\n";
    print QDP "log y\nres y $min $max\nlab NX off\n";

    # bottom window ratio
    print QDP "win 3\nlocation 0 0.067 1 0.4\n";
    print QDP "yplot on 3 7 11\n";
    print QDP "gap error\nres y\n";
    if ( $params{usetrigtime} ) {
        if ( $params{trigfrombat} ) {
            print QDP "lab X Time Since BAT Trigger (s)\n";
        } else {
            print QDP "lab X Time Since Trigger (s)\n";
        }
    } else {
        print QDP "lab X Time Since MET $params{trigtime} s (s)\n";
    }

    # scale x
    print QDP "win all\nlog x\nres x\n";

    # plot
    print QDP "cpd ${plotfile}$params{plotftype}\n";
    print QDP "plot\n";
    print QDP "quit\n";

    close QDP;

    chdir $cwd;
}

# sets file extension based on plot device
# and tries to use said device
sub getPlotDevice {

    my $stat = 0;
    if ( $params{plotftype} =~ /ps/ ) {
        $params{plotfext} = "ps";
    } elsif ( $params{plotftype} =~ /gif/ ) {
        $params{plotfext} = "gif";
    } elsif ( $params{plotftype} =~ /ppm/ ) {
        $params{plotfext} = "ppm";
    } elsif ( $params{plotftype} =~ /wd/ ) {
        $params{plotfext} = "xwd";
    } else {
        error( 1, "unsupported plot device $params{plotftype}\n" );
        return -1;
    }

    # try it!
    my $cwd = getcwd( );
    chdir $params{outdir} or die "failed to chdir to $params{outdir}\n";
   
    my $qdptest = basename( getTmpFile( "qdp" ) );
    open QDP, ">$qdptest" or die "failed to open file $qdptest\n";
    print QDP "1.0 1.0\n2.0 4.0\n";
    close QDP;
    
    my $testplot = basename( getTmpFile( "$params{plotfext}" ) ); 

    my $cmd = "qdp $qdptest 2>&1 > /dev/null";
    
    chat( 3, "testing qdp with command: $cmd\n" );
    open QDP, "|$cmd" or die "failed to run command: $cmd\n";
    print QDP "/null\n";
    print QDP "cpd ${testplot}$params{plotftype}\n";
    print QDP "plot\nquit\n";
    close QDP;

    if ( !-e $testplot ) {
        error( 1, "failed to create test plot using plot ",
               "device $params{plotftype}!\n" );
        $stat = -1;
    }
    unlink $testplot, $qdptest;
    chdir $cwd;
    return $stat;
}

#
# cleanup -
#
#       removes unneeded files, in reverse order of creation
#
sub cleanup {
    if ( $params{cleanspec} ) {
        push @clean, @cleanspec;
    }
    foreach my $cf ( reverse @clean ) {
        if ( -f $cf ) {
            debug( "cleaning un-needed file $cf\n" );
            unlink $cf;
        } elsif( -d $cf ) {
            debug( "cleaning un-needed dir $cf\n" );
            rmtree( $cf );
        }
    }
}

#
# quit -
#
sub quit {
    my $stat = shift;
    $! = $stat unless $stat == 0;
    outtro( $stat );
    #cleanup( );
    die "Task died an unnatural death\n";
}

#
# outtro
#
sub outtro {
    my $stat = shift;
    my $time  = localtime( time );
    my $success = $stat == 0 ? "success" : "failure";
    my $intro = "---- $taskName v$taskVersion $success at $time ----";
    my $dashs = '-' x length( $intro );
    chat( 1, "$dashs\n$intro\n$dashs\n\n" );
}