#! /bin/sh
# This is the LHEA perl script: /cvmfs/extras-fp7.egi.eu/extras/heasoft/heagen/x86_64-unknown-linux-gnu-libc2.19-0/bin/arftable.pl
# 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/heagen/x86_64-unknown-linux-gnu-libc2.19-0/bin/arftable.pl."
  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/heagen/x86_64-unknown-linux-gnu-libc2.19-0/bin/arftable.pl."
  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
#
# File name: arftable.pl
# Author:    Kristin Rutkowski
# $Date: 2016/03/23 01:56:11 $
# Version: 1.0
#
# This Perl script will create ARF file from photon history file
#
# Tool dependencies:
#   ftcreate
#   ftcopy
#   ftappend
#   fthedit
#
# Library dependencies:
#   gen/lib/perl/ahapp
#   gen/lib/perl/ahgen
#   gen/lib/perl/ahlog
#
# Modification history:
#
#  Ver   Date        Author  Description
#  1.0   2015-08-10  KLR     New header format. 
#



# note: offaxis angles are referred to as theta, azimuthal angles as phi


#########################
#  Pragmas
#########################

use strict;
use warnings;



#########################
#  Packages
#########################

# empty parenthesis force explicit use of namespaces
use ahlog ();
use ahgen ();
use ahapp ();           # $clobber $chatter $logfile $debug $history $mode

use File::Copy;
use List::Util qw(first);

use Data::Dumper;       # to output arrays for debugging



#########################
#  Startup
#########################

my $nargs = scalar @ARGV ;

# Query canonical APE parameters and start logging
ahapp::startup () ;



#########################
#  Input Parameters
#########################

# Query non-standard APE parameters.

# +++ 20150916 KLR should I query for vignetting before outfile, and if vignetting!=none then don't ask for outfile or angles?  Set them to blank.
my $infile           = ahapp::query_parameter("infile") ;
my $outfile          = ahapp::query_parameter("outfile") ;
my $tgtOffaxis       = ahapp::query_parameter("tgtOffaxis");
my $tgtAzimuth       = ahapp::query_parameter("tgtAzimuth");
my $vignettingfile   = ahapp::query_parameter("vignetting");

my $doVignetting     = (lc($vignettingfile) ne "none");

my $outputFile = $doVignetting ? $vignettingfile : $outfile;


#########################
#  Other Variables 
#########################

my $status = 0;

# name of extensions in history file
my $extension1 = "PHOTONHISTORY";
my $extension2 = "INPUTPHOTONS";

# name of extension in output arf file
my $extname     = "EFFECTIVEAREA";
my $vignextname = "VIGNETTING";

# history file stores up to 8 individual path codes
my $maxPathCodes          = 8;
  
# Normalization radius in mm.  This may be requested as an input param later.
# This is written to the keyword NORMRAD in the output primary HDU
my $normRadmm             = 1e30;

# precision to which to compare floating point numbers
my $precision = 16;

# array with input files
my @historyFiles;         

# how many file we will be using
my $numFilesToUse;
my @filesToUse;

# hash connecting each theta to the history file that has it
my %historyFileOfOffaxis = ();

# hash specifying which index each theta will have.  This is specific to the 
# history file with that theta, so that numInputPhotons can be accessed by 
# [fileIdx][offaxisIdx][azimIdx][energyIdx]
my %indexOfOffaxis      = ();

# total number of unique offaxis angles
my $numUnqOffaxis       = 0;

# these will be temp hashes for each file, to hold the values for that file
my %currOffaxisHash     = ();
my %currAzimHash        = ();
my %currEnergiesHash    = ();

# each history file needs to have same energy and azimuth grids
my %unqEnergieskevHash  = ();
my $numUnqEnergies      = 0;
my @currUnqEnergies     = ();
my %unqAzimDegHash      = ();
my $numUnqAzimuth       = 0;

# hashes of azimuth and energy
my %indexOfAzim         = ();
my %indexOfEnergies     = ();

# to access numInputPhotons by iFile then iTheta, need to store each theta 
# for each file.
# thetaByFile[iFile][iTheta]
my %thetaByFile;

# variables for determining EA for a normal run (not vignetting)
my $theta_nIdx          = 0;
my $theta_np1Idx        = 0;
my $azim_nIdx           = 0;
my $azim_np1Idx         = 0;
my $theta_n             = 0;
my $theta_np1           = 0;
my $azim_n              = 0;
my $azim_np1            = 0;
my $theta_nFileIdx      = 0;
my $theta_np1FileIdx    = 0;
my $theta_nIdxInFile    = 0;
my $theta_np1IdxInFile  = 0;

# array of how many input photons for each file/offaxis/azimuth/energy
# This is direct input from the second extension of each history file
my @numInputPhotons;

# each history file needs to come from same TDF, and use same geom. area
my $initTDFN            = 0;
my $initGeomAreasqcm    = 0;

# booleans checked in the main code section (doWork)
my $isDesiredAngle      = 0;
my $isTheta_n           = 0;
my $isTheta_np1         = 0;
my $isPhi_n             = 0;
my $isPhi_np1           = 0;

# output from checkPathCode()
my $isDblRefl           = 0;
my $numEightInteractions= 0;



#########################
#  Initializing 
#########################

ahlog::ah_out "Start of tool arftable.pl" ;

# Pre-Processing
$status = ahapp::begin_processing();
unless ( $status == 0 ) {
  ahlog::ah_err "BeginProcessing" ;
  ahapp::end_processing(1);
}

# Write all parameters to the log file.
ahlog::ah_info "HIGH", ahapp::write_parameters () ;

# if the user wants vignetting, only do that option.  ie: don't create 
# normal ARF.  This function will calculate that, and create the output,
# then will exit the program.
if ($doVignetting) {
} # end-if vignetting

# Check if the output file already exists.  Unless clobber is set, this will
# cause the script to fail.
unless ($ahapp::clobber) {
  if  (-e $outputFile) {
    ahlog::ah_err "Output file $outputFile already exist but clobber was not set.  Exiting." ;
    ahapp::end_processing(1);
  }
}

# If $outfile already exists, delete it.  We have already checked that clobber
# has been set.
unlink $outputFile if -e $outputFile ;

# Make sure each infile exists, before getting its data
# get list of files from infile parameter
# if the param starts with '@' (for @filelist), make sure that file exists
if ($infile =~ /^@/ and ! -e substr($infile,1)) {
  ahlog::ah_err "File $infile does not exist.";
  ahapp::end_processing(1);
}
@historyFiles = ahgen::get_file_list($infile);
foreach my $file (@historyFiles) {
  if (! -e $file) {
    ahlog::ah_err "File $file does not exist. Exiting.";
    ahapp::end_processing(1);
  }
}

my $iFile = 0;
# now that we know each input file exists, get the data
foreach my $file (@historyFiles) {

  ahlog::ah_info "HIGH", "Processing file: $file";
  
  # create a copy of the history file with ftselect
  my $tmpHistory = $file . ".tmp";
  # add to temp files, to delete later
  ahapp::add_temp_file($tmpHistory);
  
  # add the ! to clobber the temp outfile
  $status = ahgen::run_ftool("ftcopy",
                             "infile=$file",
                             "outfile=!$tmpHistory",
                             "copyall=YES",
                             "clobber=yes");
  
  # Check that photons were not originally from a file.  If they were from a 
  # file, there could have been multiple sources, and we don't want that. 
  my $srcmdl = ahgen::get_keyword($tmpHistory, $extension1, "SRCMDL");
  if ( (uc $srcmdl eq "PHOTONLIST") || (uc $srcmdl eq "GROUNDMODEL")) {
   ahlog::ah_err "History file $file must be made from photons not from a file.  When running xrtraytrace, use 'source' of point, flatcircle, or betamodel.";
    ahlog::ah_err ahgen::get_tool_stderr;
    ahapp::end_processing(1);
  }

  # Find the number of rows in second extension
  my $numRows = ahgen::get_keyword($tmpHistory, $extension2, "NAXIS2") ;

  # current number of unique offaxis theta angles, azimuthal (roll) phi angles, and energies
  my $currNumUnqOffaxis   = ahgen::get_keyword($tmpHistory, $extension2, "NOFFAXIS");
  my $currNumUnqAzimuth   = ahgen::get_keyword($tmpHistory, $extension2, "NAZIMUTH");
  my $currNumUnqEnergies  = ahgen::get_keyword($tmpHistory, $extension2, "NUMENRG");

  ahlog::ah_debug "NOFFAXIS=$currNumUnqOffaxis NAZIMUTH=$currNumUnqAzimuth NUMENRG=$currNumUnqEnergies" ;

  # grab each column into an array
  my @history2InitialTheta    = ahgen::read_column($tmpHistory, $extension2, "INITIALTHETA");
  my @history2InitialAzimDir  = ahgen::read_column($tmpHistory, $extension2, "INITIALAZIMDIR");
  my @history2Energy          = ahgen::read_column($tmpHistory, $extension2, "ENERGY");
  my @history2NumPhotons      = ahgen::read_column($tmpHistory, $extension2, "NUMPHOTONS");

  # clear the current hashes, so we can refill it with data from this file
  %currOffaxisHash    = ();
  %currAzimHash       = ();
  %currEnergiesHash   = ();

  # add the theta values to historyFileOfOffaxis hash, to record that they came 
  # from this history file (file index).  This will overwrite any duplicates, so 
  # we only have unique thetas as keys to the hash.
  for (my $iRow = 0 ; $iRow < $numRows ; $iRow++) {

    my $currOffaxis     = $history2InitialTheta[$iRow];
    my $currAzim        = $history2InitialAzimDir[$iRow];
    my $currEnergy      = $history2Energy[$iRow];
  
    $currOffaxisHash  { $currOffaxis }    = $iFile; 
    $currAzimHash     { $currAzim }       = $iFile; 
    $currEnergiesHash { $currEnergy }     = $iFile;
    
#       ahlog::ah_debug "$history2InitialTheta[$iRow] $history2InitialAzimDir[$iRow] $history2Energy[$iRow] $history2NumPhotons[$iRow] \n";

  } # end-loop through history file 2nd ext rows

  # create a sorted array of all the unique offaxis theta angles in this file
  #my @currUnqOffaxisArcmin = keys %currOffaxisHash;
  #@currUnqOffaxisArcmin= sort @currUnqOffaxisArcmin;
  my @currUnqOffaxisArcmin;
  foreach my $key (sort { $a <=> $b} keys %currOffaxisHash) {
    push @currUnqOffaxisArcmin, $key ;
  }
  
  # create a sorted array of all the unique azimuthal angles in this file
  my @currUnqAzimDeg;
  foreach my $key (sort { $a <=> $b} keys %currAzimHash) {
    push @currUnqAzimDeg, $key ;
  }
  
  # create a sorted array of all the unique energies in this file
  foreach my $key (sort { $a <=> $b} keys %currEnergiesHash) {
    push @currUnqEnergies, $key ;
  }
  
  if (ahlog::getdebug()) {
    print Data::Dumper->Dump( [ \@currUnqOffaxisArcmin ], [ qw(*currUnqOffaxisArcmin) ] );
    print Data::Dumper->Dump( [ \@currUnqAzimDeg ], [ qw(*currUnqAzimDeg) ] );
    print Data::Dumper->Dump( [ \@currUnqEnergies ], [ qw(*currUnqEnergies) ] );
  }

  # make sure that multiple files don't have the same offaxis angles
  if ($iFile > 0) {
    for (my $iOff = 0 ; $iOff < $currNumUnqOffaxis ; ++$iOff) {
      if ( exists($historyFileOfOffaxis{$currUnqOffaxisArcmin[$iOff]}) ) {
        ahlog::ah_err "A single offaxis angle must be present in one history file only.  Multiple history files cannot contain the same offaxis angle.  Offaxis angle $currUnqOffaxisArcmin[$iOff] is in multiple files";
        ahlog::ah_err ahgen::get_tool_stderr;
        ahapp::end_processing(1);
      }
    } #end-loop through unique thetas in this file
  } # end-if checking thetas in multiple files

  # now that we've checked that the current theta values aren't in 
  # historyFileOfOffaxis yet, add them
  for (my $iOff = 0 ; $iOff < $currNumUnqOffaxis ; ++$iOff) {
    $historyFileOfOffaxis{$currUnqOffaxisArcmin[$iOff]} = $iFile;
  }
  # +++ or? @historyFileOfOffaxis{ keys %currOffaxis } = @currOffaxis{ keys %currOffaxis };

  # make sure all files have the same azimuthal and energy grids, and 
  # geometric area and TDF
  if ($iFile == 0) {
    # we're at the first history file.  Store the grids.

    $numUnqEnergies = $currNumUnqEnergies;
    %unqEnergieskevHash = %currEnergiesHash;

    $numUnqAzimuth = $currNumUnqAzimuth;
    %unqAzimDegHash = %currAzimHash;

    # store the index of the azimuth and energy
    # the indices for azimuth and energy need only be stored at the first file, 
    # since every file should have the same grids.
    # these will be used for access to numInputPhotons
    for (my $iAzim = 0 ; $iAzim < $currNumUnqAzimuth ; ++$iAzim) {
      $indexOfAzim{$currUnqAzimDeg[$iAzim]} = $iAzim;
    } 
    for (my $iEnergy = 0 ; $iEnergy < $currNumUnqEnergies ; ++$iEnergy) {
      $indexOfEnergies{$currUnqEnergies[$iEnergy]} = $iEnergy;
    }
    
    # will need to make sure each history file has the same aperture geometric 
    # area and used the same TDF.  Check those keywords later as loop through 
    # history files, comparing against the keywords in the first file
    $initGeomAreasqcm  = ahgen::get_keyword($file, $extension1, "GEOMAREA");
    $initTDFN          = ahgen::get_keyword($file, $extension1, "TDFN");
    
  } else {
    # we're at the next history file.  make sure its grids match the first

    # make sure azimuth hashes have the same number of keys
    if ( keys(%currAzimHash) != keys(%unqAzimDegHash) ) {
      ahlog::ah_err "The same azimuthal grid must be used for each file.";
      ahlog::ah_err ahgen::get_tool_stderr;
      ahapp::end_processing(1);
    }
    for (my $iAzim = 0 ; $iAzim < $currNumUnqAzimuth ; ++$iAzim) {
      if ( !exists($unqAzimDegHash{$currUnqAzimDeg[$iAzim]}) ) {
        ahlog::ah_err "The same azimuthal grid must be used for each file.";
        ahlog::ah_err ahgen::get_tool_stderr;
        ahapp::end_processing(1);
      }
    } #end-loop through unique azimuthal angles in this file

    # make sure energy hashes have the same number of keys
    if ( keys(%currEnergiesHash) != keys(%unqEnergieskevHash) ) {
      ahlog::ah_err "The same energy grid must be used for each file.";
      ahlog::ah_err ahgen::get_tool_stderr;
      ahapp::end_processing(1);
    }
    for (my $iEnergy = 0 ; $iEnergy < $currNumUnqEnergies ; ++$iEnergy) {
      if ( !exists($unqEnergieskevHash{$currUnqEnergies[$iEnergy]}) ) {
        ahlog::ah_err "The same energy grid must be used for each file.";
        ahlog::ah_err ahgen::get_tool_stderr;
        ahapp::end_processing(1);
      }
    } #end-loop through unique energies in this file
    
    # +++ 20150304 KLR test these (GEOMAREA TDFN)
    #                   move this in the TRF to here
    
    # make sure every file has the same geometric area as the first one
    my $currGeomArea = ahgen::get_keyword($file, $extension1, "GEOMAREA");
    if ($currGeomArea != $initGeomAreasqcm) {
      ahlog::ah_err "All input history files must have been created using the same aperture (GEOMAREA keywords must all match).";
      ahlog::ah_err ahgen::get_tool_stderr;
      ahapp::end_processing(1);
    }
    
    # make sure every file used the same TDF as the first one we found
    my $currTDFN = ahgen::get_keyword($file, $extension1, "TDFN");
    if ($currTDFN ne $initTDFN) {
      ahlog::ah_err "All input history files must have been created using the same telescope data files (TDFN keywords must all match).";
      ahlog::ah_err ahgen::get_tool_stderr;
      ahapp::end_processing(1);
    }

  } # end-if checking grids in multiple files

  # +++ 20150303 Should this be inside the history file loop?
  # we need to store the index of the offaxis in each file, for access to 
  # numInputPhotons
  for (my $iOff = 0 ; $iOff < $currNumUnqOffaxis ; ++$iOff) {
    $indexOfOffaxis{$currUnqOffaxisArcmin[$iOff]} = $iOff;
  } 

  # store how many photons at each theta/roll/energy
  for (my $iRow = 0 ; $iRow < $numRows ; $iRow++) {

    my $currOffaxis     = $history2InitialTheta[$iRow];
    my $currAzim        = $history2InitialAzimDir[$iRow];
    my $currEnergy      = $history2Energy[$iRow];
    my $currNumPhotons  = $history2NumPhotons[$iRow];
  
    my $currOffaxisIdx  = $indexOfOffaxis{$currOffaxis};
    my $currAzimIdx     = $indexOfAzim{$currAzim};
    my $currEnergyIdx   = $indexOfEnergies{$currEnergy};
  
    # store how many input photons for this file/offaxis/azimuth/energy
    $numInputPhotons[$iFile][$currOffaxisIdx][$currAzimIdx][$currEnergyIdx] = $currNumPhotons;
    
#     $thetaByFile[$iFile][$currOffaxisIdx] = $currOffaxis;
    # also store which index this particular theta is in this particular file,
    # to access numInputPhotons later.  (duplicates will be overwritten)
    $thetaByFile{$currOffaxis} = $currOffaxisIdx;
    
#     ahlog::ah_debug "numInputPhotons[$iFile][$currOffaxis][$currAzim][$currEnergy] \t = [$iFile][$currOffaxisIdx][$currAzimIdx][$currEnergyIdx] \t = $numInputPhotons[$iFile][$currOffaxisIdx][$currAzimIdx][$currEnergyIdx] \n";
    # +++ test this by manually changing some numPhotons values to make sure matrix comes out right
  
  } # end-loop through rows

  # increment which file index we're at
  $iFile++;
  
} # end-loop through all input history files

ahlog::ah_info "LOW", "Read all the input history files, second extensions";

# create sorted array of unique offaxis theta angles
my @allUnqOffaxis;
foreach my $key (sort { $a <=> $b} keys %historyFileOfOffaxis) {
  push @allUnqOffaxis, $key ;
}
$numUnqOffaxis = keys(%historyFileOfOffaxis);

# create a hash of all the unique thetas, to store their index within allUnqOffaxis
my %indexOfAllUnqOffaxisHash;
for (my $iOff = 0 ; $iOff < $numUnqOffaxis ; ++$iOff) {
  $indexOfAllUnqOffaxisHash{$allUnqOffaxis[$iOff]} = $iOff;
}

# create sorted array of unique azimuth angles
my @allUnqAzimuth;
foreach my $key (sort { $a <=> $b} keys %unqAzimDegHash) {
  push @allUnqAzimuth, $key ;
}

# create a hash of all the unique azim, to store their index within allUnqAzimuth
my %indexOfAllUnqAzimuthHash;
for (my $iAzim = 0 ; $iAzim < $numUnqAzimuth ; ++$iAzim) {
  $indexOfAllUnqAzimuthHash{$allUnqAzimuth[$iAzim]} = $iAzim;
}

# Print info about files
ahlog::ah_info "LOW", "Geometric area squared = $initGeomAreasqcm";
ahlog::ah_info "LOW", "Energy Grid:";
for (my $iEnergy = 0 ; $iEnergy < $numUnqEnergies; ++$iEnergy) {
  ahlog::ah_info "LOW", "  $currUnqEnergies[$iEnergy]";
}
ahlog::ah_info "LOW", "Azimuthal Grid:";
for (my $iAzim = 0 ; $iAzim < $numUnqAzimuth ; ++$iAzim) {
  ahlog::ah_info "LOW", "  $allUnqAzimuth[$iAzim]";
}
ahlog::ah_info "LOW", "Offaxis Grid:";
for (my $iOff = 0 ; $iOff < $numUnqOffaxis ; ++$iOff) {
  ahlog::ah_info "LOW", "  $allUnqOffaxis[$iOff]";
}

if (ahlog::getdebug()) {
  print Data::Dumper->Dump( [ \@allUnqOffaxis ], [ qw(*allUnqOffaxis) ] );
  print Data::Dumper->Dump( [ \@allUnqAzimuth ], [ qw(*allUnqAzimuth) ] );
}

# if there was only one input angle (offaxis or azimuth), and if it does not 
# equal the target angle, the code cannot complete
if ( ($numUnqOffaxis == 1) && ($allUnqOffaxis[0] != $tgtOffaxis) ) {
    ahlog::ah_err "There is only one unique offaxis angle in the input file ($allUnqOffaxis[0]).  This must match the target offaxis angle ($tgtOffaxis).  Exiting." ;
    ahapp::end_processing(1);
}
if ( ($numUnqAzimuth == 1) && ($allUnqAzimuth[0] != $tgtAzimuth) ) {
    ahlog::ah_err "There is only one unique azimuthal angle in the input file ($allUnqAzimuth[0]).  This must match the target offaxis angle ($tgtAzimuth).  Exiting." ;
    ahapp::end_processing(1);
}

# if we're doing vignetting, then we need all the input files
if ($doVignetting) {

  @filesToUse = @historyFiles;

} else {

# Determine the four surrounding theta/phi values, where theta_0 and phi_0 are 
  # the target offaxis and roll, and n is the index for the preceding values, 
  # and n+1 is the index for the values after. 
  # These are the indices we'll need for the surrounding offaxis and azimuth 
  # angles.  The offaxis is respective to the history file in which it's found.
  # These will be needed for access to numInputPhotons
  #     if phi_0 == target, then phi_n = phi_0 = phi_n+1 (same for theta)
  #     phi_n <= phi_0 <= phi_n+1
  #     theta_n <= theta_0 <= theta_n+1
  ($theta_nIdx, $theta_np1Idx) = bisectionLocate(\@allUnqOffaxis, $numUnqOffaxis, $tgtOffaxis);
  ($azim_nIdx, $azim_np1Idx) = bisectionLocate(\@allUnqAzimuth, $numUnqAzimuth, $tgtAzimuth);

  # NOTE: If the target offaxis (or azimuth) is exactly on the grid of input 
  # offaxis, we don't technically need to call bisectionLocate or interpolate. 
  # Doing it this way will make ea1 and ea2 (for example) equal, and make 
  # interpolate inefficient. But it will make code ready to distribute quicker.

  ahlog::ah_debug "theta_nIdx=$theta_nIdx theta_np1Idx=$theta_np1Idx \n";
  ahlog::ah_debug "azim_nIdx=$azim_nIdx azim_np1Idx=$azim_np1Idx \n";

  # store the offaxis and azimuthal angles we will need
  $theta_n     = $allUnqOffaxis[$theta_nIdx]; 
  $theta_np1   = $allUnqOffaxis[$theta_np1Idx];
  $azim_n      = $allUnqAzimuth[$azim_nIdx];
  $azim_np1    = $allUnqAzimuth[$azim_np1Idx];

  # note which file(s) we will need to read based on the target offaxis and roll.
  # even if both thetas are in the same file, this will put that same filename
  # into both entries
  $theta_nFileIdx = $historyFileOfOffaxis{$theta_n};
  $theta_np1FileIdx = $historyFileOfOffaxis{$theta_np1};
  $filesToUse[0] = $historyFiles[$theta_nFileIdx];
  if ($theta_nFileIdx != $theta_np1FileIdx) {
    $filesToUse[1] = $historyFiles[$theta_np1FileIdx];
  }

  # +++ new 20150325.  add to trf
  # now that we know the theta indices, and which file each comes from, we need
  # to get the theta index particular to that file, so that numInputPhotons
  # can be accessed by file index then by theta index
  $theta_nIdxInFile    = $thetaByFile{$theta_n};
  $theta_np1IdxInFile  = $thetaByFile{$theta_np1};

  ahlog::ah_debug "theta_n   = %4f (index %d) (theta_nFileIdx=%d, theta_nIdxInFile=%d) \n",     $theta_n,   $theta_nIdx,    $theta_nFileIdx,    $theta_nIdxInFile;
  ahlog::ah_debug "theta_np1 = %4f (index %d) (theta_np1FileIdx=%d, theta_np1IdxInFile=%d) \n",  $theta_np1, $theta_np1Idx,  $theta_np1FileIdx,  $theta_np1IdxInFile  ;
  ahlog::ah_debug "azim_n    = %4f (index %d) \n", $azim_n, $azim_nIdx;
  ahlog::ah_debug "azim_np1  = %4f (index %d) \n", $azim_np1, $azim_np1Idx;

} # end-if vignetting

$numFilesToUse = scalar @filesToUse;

ahlog::ah_info "HIGH", "Number of files to use = $numFilesToUse";
ahlog::ah_info "HIGH", "Using history file(s): @filesToUse";

ahlog::ah_debug "filesToUse: @filesToUse \n";

# grab other keywords from the first history file, to write to output later
my $telescop          = ahgen::get_keyword($filesToUse[0], $extension1, "TELESCOP");
my $instrume          = ahgen::get_keyword($filesToUse[0], $extension1, "INSTRUME");
my $detnam            = ahgen::get_keyword($filesToUse[0], $extension1, "DETNAM");
my $cvsd0001          = ahgen::get_keyword($filesToUse[0], $extension1, "CVSD0001");
my $cvst0001          = ahgen::get_keyword($filesToUse[0], $extension1, "CVST0001");

# set up effective area (EA) arrays for the four surrounding theta/phi 
# combinations, the interpolated EA, and the final EA.  initialize to 0
# EA_1(E) for theta_n, phi_n
# EA_2(E) for theta_n, phi_n+1
# EA_3(E) for theta_n+1, phi_n
# EA_4(E) for theta_n+1, phi_n+1
my @ea1                 = (0) x $numUnqEnergies;
my @ea2                 = (0) x $numUnqEnergies; 
my @ea3                 = (0) x $numUnqEnergies; 
my @ea4                 = (0) x $numUnqEnergies;
my @ea_phi_0_theta_n    = (0) x $numUnqEnergies;
my @ea_phi_0_theta_np1  = (0) x $numUnqEnergies;
my @eaFinal             = (0) x $numUnqEnergies;

# these are for vignetting, if that is requested
# +++ 20150804 KLR change of plans - we're not averaging over azimuth - add to TRF
my @eaVignette          = ();

# This will hold the final vignette ea ratio.  We can't write over the 
# original, becuase then the values for theta=0 will be written as 1 after 
# it divides by itself.
my @eaVignetteFinal = ();

# now that we know number of unique offaxis, azimuth, energy: initialize arrays
# (numOffaxis, numAzim, numUnqEnergies)
for (my $iTheta = 0 ; $iTheta < $numUnqOffaxis; ++$iTheta) {
  for (my $iEnergy = 0 ; $iEnergy < $numUnqEnergies; ++$iEnergy) {
    for (my $iAzim = 0 ; $iAzim < $numUnqAzimuth; ++$iAzim) {
      $eaVignette[$iTheta][$iAzim][$iEnergy] = 0;
    }
    #$eaVigSumAzim[$iTheta][$iEnergy] = 0;
  }
}

# now that we've initizalized
my $foundZero = 0;
if ($doVignetting) {
  ahlog::ah_info "LOW", "Checking that theta=0 is one of the inputs";
  # make sure that theta=0 is one of the inputs
  for (my $iTheta = 0 ; $iTheta < $numUnqOffaxis; ++$iTheta) {
    if (equal($allUnqOffaxis[$iTheta], 0, $precision)) {
      $foundZero = 1;
      last;
    }
  }
  if (!$foundZero) {
    ahlog::ah_err "Offaxis=0 must be in one of the input history files in order to create a vignetting file.";
    ahlog::ah_err ahgen::get_tool_stderr;
    ahapp::end_processing(1);
    
  }
}


#########################
#  Main Code Block 
#########################

# read files we decided to keep based on which file has which theta and roll
foreach my $file (@filesToUse) {

  ahlog::ah_info "HIGH", "Reading input file: $file";

  # create a copy of the history file with ftselect
  my $tmpHistory = $file . ".main.tmp";
  # add to temp files, to delete later
  ahapp::add_temp_file($tmpHistory);
  
  # add the ! to clobber the temp outfile
  $status = ahgen::run_ftool("ftcopy",
                             "infile=$file",
                             "outfile=!$tmpHistory",
                             "copyall=YES",
                             "clobber=yes");
  
  # find out if history file contains only photons that hit the results plane
  my $resultsPlaneOnly = ahgen::get_keyword($tmpHistory, $extension1, "RSLPONLY");

  # if the history file contained only photons that hit the results plane, 
  # then we don't need to check for the resultsplane column.
  my $checkForResultsPlane = !$resultsPlaneOnly;

  # store how many total photons (rows) in current history file
  my $numPhotons = ahgen::get_keyword($tmpHistory, $extension1, "NAXIS2");

  # store the following history file columns in local arrays
  my @historyEnergy               = ahgen::read_column($tmpHistory, $extension1, "energy");
  my @historyInitialThetaArcmin   = ahgen::read_column($tmpHistory, $extension1, "initialtheta");
  my @historyInitialAzimDirDeg    = ahgen::read_column($tmpHistory, $extension1, "initialazimdir");
  my @historyPathcode             = ahgen::read_column($tmpHistory, $extension1, "pathcode");
  my @historyResultsPlane;
  if ($checkForResultsPlane) {
    @historyResultsPlane          = ahgen::read_column($tmpHistory, $extension1, "resultsplane");
  } else {
    # else, everything in file did hit results plane. initialize all true
    @historyResultsPlane = ((1) x $numPhotons);
  }

  # loop through rows of current history file
  for (my $iRow = 0 ; $iRow < $numPhotons ; $iRow++) {
  
    ahlog::ah_debug " **** iRow:$iRow theta:$historyInitialThetaArcmin[$iRow] roll:$historyInitialAzimDirDeg[$iRow] pathcode:$historyPathcode[$iRow] resultsplane:$historyResultsPlane[$iRow] \n";

    # reset these booleans to 0 before we go into checkForAzimOffaxis
    $isTheta_n = 0;
    $isTheta_np1 = 0;
    $isPhi_n = 0;
    $isPhi_np1 = 0;
    
    # if we're on a normal run (not vignetting), check for which angle we're at
    if (!$doVignetting) {
      # only keep history file rows with the desired phi and theta values
      ($isDesiredAngle, $isTheta_n, $isTheta_np1, $isPhi_n, $isPhi_np1) = 
          checkForAzimOffaxis($historyInitialThetaArcmin[$iRow], $historyInitialAzimDirDeg[$iRow]);
      if (!$isDesiredAngle) {
        next;
      }
    }
    
    # call the function to see if the current photon hit the results plane and 
    # was a double reflection.  If it wasn't, move on to the next photon.
    ($isDblRefl, $numEightInteractions) = checkPathCode($historyPathcode[$iRow], $historyResultsPlane[$iRow]);
    if (!$isDblRefl) {
      next;
    }
    
    # if we've made it this far, we know this photon should be counted towards 
    # the appropriate effective area
    
    # figure out which energy index we're at
    my $currEnergy = $historyEnergy[$iRow];
    my $currEnergyIdx = $indexOfEnergies{$currEnergy};  
    #ahlog::ah_debug "currEnergyIdx of $currEnergy = $currEnergyIdx \n";

    # if we're doing vignetting, we need to use all the angles
    if ($doVignetting) {
      
      # figure out which theta and azim index we're at
      my $currOff = $historyInitialThetaArcmin[$iRow];
      #my $currOffIdx = $indexOfOffaxis{$currOff};
      my $currOffIdx = $indexOfAllUnqOffaxisHash{$currOff};
      
      my $currAzim = $historyInitialAzimDirDeg[$iRow];
      #my $currAzimIdx = $indexOfAzim{$currAzim};
      my $currAzimIdx = $indexOfAllUnqAzimuthHash{$currAzim};
    
      ahlog::ah_debug "calculating vignetting currOff=$currOff currOffIdx=$currOffIdx currAzim=$currAzim currAzimIdx=$currAzimIdx \n";
      
      $eaVignette[$currOffIdx][$currAzimIdx][$currEnergyIdx]++;
      
#      ahlog::ah_debug "eaVignette[$iTheta][$iAzim][$iEnergy] = $eaVignette[$iTheta][$iAzim][$iEnergy] \n";
      
    } else {
      # if we're on a normal run (not vignetting), calculate approp EAs
    
      ahlog::ah_debug "calculating normal ea \n";
      
      # now fill in effective area (EA) arrays for the four surrounding theta/phi 
      # combinations, whichever theta and phi this current row is
      # ea1(E) for theta_n, phi_n
      # ea2(E) for theta_n, phi_n+1
      # ea3(E) for theta_n+1, phi_n
      # ea4(E) for theta_n+1, phi_n+1

      # figure out which ea array this theta/phi combo belongs in
      # +++ 20150804 KLR This isn't the most efficient, if a target angle is exactly on the grid.
      #                  I probably should avoid interpolating at all if that is the case
      if ($isTheta_n) {
        if ($isPhi_n) {
          $ea1[$currEnergyIdx]++;
          #ahlog::ah_debug "isTheta_n and isPhi_n: ea1[$currEnergyIdx] = $ea1[$currEnergyIdx] \n";
        } 
        if ($isPhi_np1) {
          $ea2[$currEnergyIdx]++;
          #ahlog::ah_debug "isTheta_n and isPhi_np1: ea2[$currEnergyIdx] = $ea2[$currEnergyIdx] \n";
        }
      } 
      if ($isTheta_np1) {
        if ($isPhi_n) {
          $ea3[$currEnergyIdx]++;
          #ahlog::ah_debug "isTheta_np1 and isPhi_n: ea3[$currEnergyIdx] = $ea3[$currEnergyIdx] \n";
        }
        if ($isPhi_np1) {
          $ea4[$currEnergyIdx]++;
          #ahlog::ah_debug "isTheta_np1 and isPhi_np1: ea4[$currEnergyIdx] = $ea4[$currEnergyIdx] \n";
        }
      }
      
    } # end-if we're doing vignetting or not
  
  } # end-loop through rows in history file

} #end-loop through history files

ahlog::ah_info "HIGH", "Finished reading all input files.";

# alert the user if there were any pathcodes with 8 interactions.  Since the 
# history file only records 8 interactions, this photon may have interactions 
# we're missing.
if ($numEightInteractions > 1) {
  ahlog::ah_info "LOW", "There were $numEightInteractions photons with at least $maxPathCodes interactions.";
} elsif ($numEightInteractions > 0) {
  # there was only one.
  ahlog::ah_info "LOW", "There was $numEightInteractions photon with at least $maxPathCodes interactions.";
}

# finish calculating the EA now that we’ve counted the photons
# EA = (numimpactphotons/numtotalphotons) * geometric area
    
# if we're doing vignetting, we need to use all the angles
if ($doVignetting) {

  ahlog::ah_info "HIGH", "Calculating EA for vignetting file";

  for my $iTheta (0 .. $#allUnqOffaxis) {
    my $currTheta = $allUnqOffaxis[$iTheta];
    # get the indices needed for this current theta
    my $theta_FileIdx     = $historyFileOfOffaxis{$currTheta};
    my $theta_IdxInFile   = $thetaByFile{$currTheta};

    for (my $iAzim = 0 ; $iAzim < $numUnqAzimuth ; ++$iAzim) {
      for (my $iEnergy = 0 ; $iEnergy < $numUnqEnergies ; ++$iEnergy) {
        $eaVignette[$iTheta][$iAzim][$iEnergy] /= $numInputPhotons[$theta_FileIdx][$theta_IdxInFile][$iAzim][$iEnergy];
        $eaVignette[$iTheta][$iAzim][$iEnergy] *= $initGeomAreasqcm;
      } # end-loop through energy
    } # end-loop through azimuth
  } # end-loop through offaxis

  my $theta0Idx = first { $allUnqOffaxis[$_] == 0.0 } 0..$#allUnqOffaxis;
  ahlog::ah_debug "theta0Idx=$theta0Idx \n";
  
  # now we can divide by EA_0 (EA at theta=0)
  for my $iTheta (0 .. $#allUnqOffaxis) {
    for (my $iEnergy = 0 ; $iEnergy < $numUnqEnergies ; ++$iEnergy) {
      for (my $iAzim = 0 ; $iAzim < $numUnqAzimuth ; ++$iAzim) {
        ahlog::ah_debug "eaVignette[$iTheta][$iAzim][$iEnergy] / eaVignette[$theta0Idx][$iAzim][$iEnergy] = $eaVignette[$iTheta][$iAzim][$iEnergy]/$eaVignette[$theta0Idx][$iAzim][$iEnergy]  \n";
        $eaVignetteFinal[$iTheta][$iAzim][$iEnergy] = $eaVignette[$iTheta][$iAzim][$iEnergy] / $eaVignette[$theta0Idx][$iAzim][$iEnergy];
      }
    }
  }
  
#  # now we can sum over azimuth
#  for my $iTheta (0 .. $#allUnqOffaxis) {
#    for (my $iEnergy = 0 ; $iEnergy < $numUnqEnergies ; ++$iEnergy) {
#      for (my $iAzim = 0 ; $iAzim < $numUnqAzimuth ; ++$iAzim) {
#        $eaVigSumAzim[$iTheta][$iEnergy] += $eaVigByAzim[$iTheta][$iAzim][$iEnergy];
#      }
#    }
#  }
#
#  # now we can divide by EA_0 (EA at theta=0)
#  for my $iTheta (0 .. $#allUnqOffaxis) {
#    for (my $iEnergy = 0 ; $iEnergy < $numUnqEnergies ; ++$iEnergy) {
#      $eaVigSumAzim[$iTheta][$iEnergy] /= $eaVigSumAzim[$theta0Idx][$iEnergy];
#    }
#  }
  
  if (ahlog::getdebug()) {
    ahlog::ah_debug "After all EA Vignetting calculations \n";
    print Data::Dumper->Dump( [ \@eaVignette ], [ qw(*eaVignette) ] );
  }

} else {
  # we're not doing vignetting.  
  
  ahlog::ah_info "HIGH", "Calculating EA for target offaxis and target azimuth";

  if (ahlog::getdebug()) {
    print Data::Dumper->Dump( [ \@ea1 ], [ qw(*ea1) ] );
    print Data::Dumper->Dump( [ \@ea2 ], [ qw(*ea2) ] );
    print Data::Dumper->Dump( [ \@ea3 ], [ qw(*ea3) ] );
    print Data::Dumper->Dump( [ \@ea4 ], [ qw(*ea4) ] );
  }

  for (my $iEnergy = 0 ; $iEnergy < $numUnqEnergies ; ++$iEnergy) {

    #   ahlog::ah_debug " ***** iEnergy $iEnergy of $numUnqEnergies \n";
    #   ahlog::ah_debug "ea1[$iEnergy] = $ea1[$iEnergy] \n";
    #   ahlog::ah_debug "ea2[$iEnergy] = $ea2[$iEnergy] \n";
    #   ahlog::ah_debug "ea3[$iEnergy] = $ea3[$iEnergy] \n";
    #   ahlog::ah_debug "ea4[$iEnergy] = $ea4[$iEnergy] \n";

    #   ahlog::ah_debug "numInputPhotons[$theta_nFileIdx][$theta_nIdxInFile][$azim_nIdx][$iEnergy] = $numInputPhotons[$theta_nFileIdx][$theta_nIdxInFile][$azim_nIdx][$iEnergy] \n";
    #   ahlog::ah_debug "numInputPhotons[$theta_nFileIdx][$theta_nIdxInFile][$azim_np1Idx][$iEnergy] = $numInputPhotons[$theta_nFileIdx][$theta_nIdxInFile][$azim_np1Idx][$iEnergy] \n";
    #   ahlog::ah_debug "numInputPhotons[$theta_np1FileIdx][$theta_np1IdxInFile][$azim_nIdx][$iEnergy] = $numInputPhotons[$theta_np1FileIdx][$theta_np1IdxInFile][$azim_nIdx][$iEnergy] \n";
    #   ahlog::ah_debug "numInputPhotons[$theta_np1FileIdx][$theta_np1IdxInFile][$azim_np1Idx][$iEnergy] = $numInputPhotons[$theta_np1FileIdx][$theta_np1IdxInFile][$azim_np1Idx][$iEnergy] \n";

    $ea1[$iEnergy] /= $numInputPhotons[$theta_nFileIdx][$theta_nIdxInFile][$azim_nIdx][$iEnergy];
    $ea2[$iEnergy] /= $numInputPhotons[$theta_nFileIdx][$theta_nIdxInFile][$azim_np1Idx][$iEnergy];
    $ea3[$iEnergy] /= $numInputPhotons[$theta_np1FileIdx][$theta_np1IdxInFile][$azim_nIdx][$iEnergy];
    $ea4[$iEnergy] /= $numInputPhotons[$theta_np1FileIdx][$theta_np1IdxInFile][$azim_np1Idx][$iEnergy];
    $ea1[$iEnergy] *= $initGeomAreasqcm;
    $ea2[$iEnergy] *= $initGeomAreasqcm;
    $ea3[$iEnergy] *= $initGeomAreasqcm;
    $ea4[$iEnergy] *= $initGeomAreasqcm;

  } #end-loop through energies in EA

  ahlog::ah_debug "finished calculating ea1 ea2 ea3 ea4 \n";
  
  if (ahlog::getdebug()) {
     print Data::Dumper->Dump( [ \@ea1 ], [ qw(*ea1) ] );
     print Data::Dumper->Dump( [ \@ea2 ], [ qw(*ea2) ] );
     print Data::Dumper->Dump( [ \@ea3 ], [ qw(*ea3) ] );
     print Data::Dumper->Dump( [ \@ea4 ], [ qw(*ea4) ] );
  }

  # now start interpolating effective areas to get the final EA

  # First: interpolate on phi (roll) at theta_n
  # ea_phi_0_theta_n(E) = interpolate btwn EA_1(E) and EA_2(E)
  @ea_phi_0_theta_n = @{interpolate(\@ea1, \@ea2, $azim_n, $azim_np1, $tgtAzimuth)};
  
  # Second: interpolate on phi (roll) at theta_n+1
  # ea_phi_0_theta_np1(E) = interpolate btwn EA_3(E) and EA_4(E)
  @ea_phi_0_theta_np1 = @{interpolate(\@ea3, \@ea4, $azim_n, $azim_np1, $tgtAzimuth)};

  # Third: interpolate on theta (off axis) to get final
  # eaFinal(E) = interpolate btwn EA_phi_0_theta_n(E) and EA_phi_0_theta_np1(E)
  @eaFinal = @{interpolate(\@ea_phi_0_theta_n, \@ea_phi_0_theta_np1, $theta_n, $theta_np1, $tgtOffaxis)};

  if (ahlog::getdebug()) {
    print Data::Dumper->Dump( [ \@ea_phi_0_theta_n ], [ qw(*ea_phi_0_theta_n) ] );
    print Data::Dumper->Dump( [ \@ea_phi_0_theta_np1 ], [ qw(*ea_phi_0_theta_np1) ] );
    print Data::Dumper->Dump( [ \@eaFinal ], [ qw(*eaFinal) ] );
    print Data::Dumper->Dump( [ \@currUnqEnergies ], [ qw(*currUnqEnergies) ] );
  }

}  # end-if we're doing vignetting or not
  
ahlog::ah_info "HIGH", "Finished Calculating EA";

# now write the output file

if ($doVignetting) {
  
  # for vignetting, we want an extension for every theta and azimuth
  
  # in Perl, using ftcreate, can only create one file/extension at a time.  Add
  # those names to this array.  Then we'll merge them all into one file using 
  # ftmerge after the loops
  my @outfiles;
  
  # first, create a file with each extension
  my $extNum = 0;
  my $extname;
  my $currOffaxis = 0.0;
  my $currAzimuth = 0.0;
  for my $iTheta (0 .. $#allUnqOffaxis) {
    $currOffaxis = $allUnqOffaxis[$iTheta];
    
    for (my $iAzim = 0 ; $iAzim < $numUnqAzimuth ; ++$iAzim) {
      $currAzimuth = $allUnqAzimuth[$iAzim];
    
      ++$extNum;
      $extname = "VIGNETTE_EA_" . $extNum;
      ahlog::ah_debug "extname = $extname \n";
      
      # current name for this file (one for each extension, until combine)
      my $currOutfile = $outputFile . ".extension" . $extNum;
      push @outfiles, $currOutfile;
      ahapp::add_temp_file($currOutfile);
      
      # put the current ea (for this theta and roll) into its own array
      my @currEA;
      for (my $iEnergy = 0 ; $iEnergy < $numUnqEnergies ; ++$iEnergy) {
        $currEA[$iEnergy] = $eaVignetteFinal[$iTheta][$iAzim][$iEnergy];
      }
      
      # write the output bintable extension
      # +++ 20150805 KLR need to pass in the array with this EA
      writeExtension($currOutfile, $extname, $currOffaxis, $currAzimuth, \@currEA);

    } # end-loop through azimuth
  } # end-loop through offaxis
  
  # grab the first file we made, and make it the regular outfile
  $status = ahgen::run_ftool("ftcopy",
                  "infile=$outfiles[0]",
                  "outfile=$outputFile",
                  "clobber=yes");
  if ($status) {
    ahlog::ah_err "Error creating file $outputFile with first extension";
    ahapp::end_processing(1);
  }
  
  # now merge all the rest of the files into the outfile
  for (my $iFile = 1 ; $iFile < $extNum ; ++$iFile) {
      $status = ahgen::run_ftool("ftappend",
                    "infile=$outfiles[$iFile]",
                    "outfile=$outputFile");
    if ($status) {
      ahlog::ah_err "Error adding extension $iFile to $outputFile";
      ahapp::end_processing(1);
    }
  }
  
} else {

  # in normal mode (not vignetting), we create a single extension for the
  # target offaxis and target azimuth

  # write the output EA bintable extension
  writeExtension($outputFile, $extname, $tgtOffaxis, $tgtAzimuth, \@eaFinal);
  
} # end-if vignetting or not



#########################
#  Finishing
#########################


ahlog::ah_debug "outputFile = $outputFile \n";

# Add CHECKSUM and DATASUM, and check that file we just created is valid
unless (ahgen::update_checksum_and_verify($outputFile)) {
  ahgen::ah_err "FITS file $outputFile failed FITS verification test." ;
  ahapp::end_processing(1);
}
ahlog::ah_info "HIGH", "Successfully created FITS file: $outputFile" ;    

# We're done
ahapp::end_processing($status);





#########################
#  Subroutines
#########################


#############################################################################

sub bisectionLocate {

  ahlog::ah_debug "---- start bisectionLocate ---- \n";

  # -------------------------------------
  # define and initialize all variables
  # -------------------------------------
  
  # input parameters
  my $xGridRef    = shift;
  my @xGrid       = @{$xGridRef};
  my $numGridPts  = shift;
  my $xIn         = shift;
  
  ahlog::ah_debug "numGridPts = $numGridPts \n";
  ahlog::ah_debug "xIn = $xIn \n";

  # initial values of pointers, they will change as the routine progresses
  my $lobin = 0;
  my $hibin = $numGridPts - 1;
  my $deltabins = $hibin - $lobin + 1;
  my $midbin = 0;
  my $foundBins = 0;
  
  my $firstX = $xGrid[$lobin];
  my $lastX  = $xGrid[$hibin];
  my $midX   = 0;
  
  ahlog::ah_debug "lobin=$lobin hibin=$hibin deltabins=$deltabins midbin=$midbin firstX=$firstX lastX=$lastX midX=$midX \n";
  
  # output parameters
  my $index1 = 0;
  my $index2 = 0;
  
  # -------------------------------------
 
 
  # case of XIN = highest grid value
  if ($xIn == $lastX) {
    $index1 = $numGridPts - 1;
    $index2 = $index1;
    return ($index1, $index2);
  }
  
  # case of XIN = lowest grid value
  if ($xIn == $firstX) {
    $index1 = 0;
    $index2 = $index1;
    return ($index1, $index2);
  }
 
  # make sure xin is valid
  if ( ($xIn < $firstX) || ($xIn > $lastX) ) {
    ahlog::ah_err "Input x=$xIn is not in range ($firstX,$lastX)";
    ahlog::ah_err ahgen::get_tool_stderr;
    ahapp::end_processing(1);
  }
  
  my $iter = 0;
  while ($deltabins > 2) {
    ahlog::ah_debug "iter=$iter: lobin=$lobin hibin=$hibin deltabins=$deltabins midbin=$midbin \n";
    
    $midbin = $lobin + int($deltabins/2);
    $midX = $xGrid[$midbin];
    if ( ($xIn > $midX) && ($xIn < $lastX) ) {
      $lobin = $midbin;
    } elsif ( ($xIn < $midX) && ($xIn > $firstX) ) {
      $hibin = $midbin;
    # Case of xIn exactly equal to the mid value
    } elsif ($xIn == $midX) {
      $index1 = $midbin - 1;
      $index2 = $midbin;
      $foundBins = 1;
      ahlog::ah_debug "foundBins\n";
      last;
    }
    # lobin and hibin move - following is the latest deltabins
    $deltabins = $hibin - $lobin + 1;

    $iter++;
  } # end-while
  
  # if we didn't find the correct bins earlier, set them to lobin and hibin
  if (!$foundBins) {
    ahlog::ah_debug "index1=lobin\n";
    # lobin and hibin should now point to the index1 and index2 respectively 
    $index1 = $lobin;
    $index2 = $hibin;
  }

  ahlog::ah_debug "---- end bisectionLocate ---- \n";
  
  return ($index1, $index2);

} # end-sub bisectionLocate


#############################################################################


# \brief Analyzes the pathcode for a single photon to determine the properties
# \param[in] pathcode The string containing the entire pathcode.  It is an 
#                 array of three integers per interaction that record various 
#                 attributes of the interaction. The three integers are 
#                 referred to as A, BC, and D, and are defined as follows:
#      pathcode[0] = A = interaction type
#        1 = absorption or end of path on the results plane
#        2 = reflection only
#        3 = reflection followed by scattering
#        4 = transmission
#        9 = error
#      pathcode[1] = BC = impacted object
#        01 = results plane
#        02 = inner housing
#        03 = outer housing
#        04 = segment/sector sidewall
#        05 = pre-collimator foil
#        06 = primary mirror foil
#        07 = secondary mirror foil
#        08 = support structure
#        09 = top external object (above main telescope)
#        10 = bottom external object (below main telescope)
#      pathcode[2] = D = impacted face
#         1 = back (mirrors, pre-collimator foils)
#             i.e. the face oriented towards outer housing
#             for obstructions & focal plane, it is the top
#         2 = front (mirrors, pre-collimator foils)
#             i.e. the face oriented towards inner housing
#             for obstructions it is the bottom.
#         3 = mirror, pre-collimator foils: top face/edge 
#         4 = mirror, pre-collimator foils: bottom face/edge
#         5 = mirror, pre-collimator foils: sides
#         6 = Undetermined (but not back or front)
# \param[in] resultsPlane True if this photon hit the results plane
# \param[out] isDblRefl True if this photon was a double reflection, meaning it 
#                 reflected off the primary mirrors exactly once and the 
#                 secondary mirrors exactly once, then hit the results plane.  
#                 These reflections may include scattering.  There may also be 
#                 transmission events in the pathcode, and events involving the 
#                 precollimator or external objects: these do not effect the 
#                 double reflection status.
# \param[out] numEightInteractions The history file only writes up to 8
#                 interactions in the pathcode.  If there were 8 portions of 
#                 the current pathcode, we don't know if there were exactly 8 
#                 or more than 8.  We will still count this photon (as long as 
#                 it meets other criteria), but we will increment a counter of 
#                 these photons and alert the user to it after looping through 
#                 all the photons.
#
# checkPathCode analyzes the pathcode of a single photon to determine if this
# was a double reflection, and if this pathcode records 8 interactions.
# 
sub checkPathCode {

#  ahlog::ah_debug "---- start checkPathCode ---- \n";

  # -------------------------------------
  # define and initialize all variables
  # -------------------------------------
  
  # input parameters
  my $pathcode              = shift;
  my $resultsPlane          = shift;

  # set some constants
  my $indivPathCodeSize = 4;      # size of an individual path code

  # strings that represent reflection pathcodes to search for later.
  # these strings are hard coded into xrtraytrace
  my $priRefl       = "2062";
  my $priReflScat   = "3062";
  my $secRefl       = "2072";
  my $secReflScat   = "3072";

  # count how many reflections from primary and secondary mirrors
  my $numPrimRefl   = 0; 
  my $numSecRefl    = 0;
  
  # output parameters
  my $isDblRefl             = 0;
  my $hasEightInteractions  = 0;

  # -------------------------------------
  
  # first, check if photon hit the results plane
  # if we need to check the resultsplane column, and it was false (meaning 
  # the photon in this row did not hit the results plane), we can return from 
  # this function and move to next photon
  if (!$resultsPlane) {
    return ($isDblRefl, $hasEightInteractions);
  }
  
  # next, check if this was a double reflection

  # Check pathcode column.  Only process double reflections, which could 
  # include transmission, scattering, and external objects.  We could accept 
  # all photons, which would be closer to real data, as real data doesn't know 
  # which was a double reflection, but we're deciding to only process double 
  # reflections.  Don't count if more than one primary + one secondary refl.

  # see if there were 8 interactions in this pathcode.  If so, we will alert
  # user later
  my $numPathCodeChars = length($pathcode);
  my $numInteractions = $numPathCodeChars / $indivPathCodeSize;
  if ($numInteractions = $maxPathCodes) {
    $hasEightInteractions++;
  }

  # break up the pathcode into individual pieces of 4
  my @currPathCodes = unpack("(A4)*", $pathcode);

  # loop through the pathcodes, looking for a double reflection
  foreach my $code (@currPathCodes) {
  
    # we only care about reflections from primary and secondary.  We can ignore
    # all other interactions. (transmission, external objects, etc)

    # search for reflections from primaries (can have scattering)
    if ( ($code == $priRefl) || ($code == $priReflScat) ) {
      $numPrimRefl++;
    }

    # seach for reflections from secondaries (can have scattering)
    if ( ($code == $secRefl) || ($code == $secReflScat) ) {
      $numSecRefl++;
    }

  } # end-loop through pathcode

  # Make sure there was exactly one primary reflection and one secondary 
  # reflection. If not, move on to the next photon.
  if ( ($numPrimRefl == 1) && ($numSecRefl == 1) ) {
    $isDblRefl = 1;
  }

#  ahlog::ah_debug "---- end checkPathCode ---- \n";
  
  return ($isDblRefl, $hasEightInteractions);
  
} # end-sub checkPathCode


#############################################################################


sub checkForAzimOffaxis {

#  ahlog::ah_debug "---- start checkForAzimOffaxis ---- \n";

  # -------------------------------------
  # define and initialize all variables
  # -------------------------------------
  
  # input parameters
  my $currOffaxis   = shift;
  my $currAzim      = shift;
  # note: theta_n, theta_np1, phi_n, phi_np1 are declared earlier using 'our'
  
  # default outputs, set to false
  my $isDesiredAngle = 0;
  my $isTheta_n = 0;
  my $isTheta_np1 = 0;
  my $isPhi_n = 0;
  my $isPhi_np1 = 0;

  # -------------------------------------
 
  if ($currOffaxis == $theta_n) {
    $isTheta_n = 1;
    $isDesiredAngle = 1;
  } 
  if ($currOffaxis == $theta_np1) {
    $isTheta_np1 = 1;
    $isDesiredAngle = 1;
  }

  if ($currAzim == $azim_n) {
    $isPhi_n = 1;
    $isDesiredAngle = 1;
  } 
  if ($currAzim == $azim_np1) {
    $isPhi_np1 = 1;
    $isDesiredAngle = 1;
  }

#  ahlog::ah_debug "---- end checkForAzimOffaxis ---- \n";
  ahlog::ah_debug "isDesiredAngle=$isDesiredAngle isTheta_n=$isTheta_n isTheta_np1=$isTheta_np1 isPhi_n=$isPhi_n isPhi_np1=$isPhi_np1 \n";
  return ($isDesiredAngle, $isTheta_n, $isTheta_np1, $isPhi_n, $isPhi_np1);

} # end-sub checkForAzimOffaxis


#############################################################################


# pass in two arrays, and the variables
sub interpolate {

  ahlog::ah_debug "---- start interpolate ---- \n";

  # -------------------------------------
  # define and initialize all variables
  # -------------------------------------
  
  # input parameters
  my $inArray1Ref    = shift;             # changing 'y' array: first
  my @inArray1       = @{$inArray1Ref};
  my $inArray2Ref    = shift;             # changing 'y' array: second
  my @inArray2       = @{$inArray2Ref};
  my $var1           = shift;             # changing 'x' value: first
  my $var2           = shift;             # changing 'x' value: second
  my $varNaught      = shift;             # changing 'x' value: target
  
  # output array
  my @outArray       = (0) x $numUnqEnergies;
  
  # -------------------------------------
  
  if (ahlog::getdebug()) {
    print Data::Dumper->Dump( [ \@inArray1 ], [ qw(*inArray1) ] );
    print Data::Dumper->Dump( [ \@inArray2 ], [ qw(*inArray2) ] );
    ahlog::ah_debug "var1=$var1 \n";
    ahlog::ah_debug "var2=$var2 \n";
  }
  
  for (my $iEnergy = 0 ; $iEnergy < $numUnqEnergies ; ++$iEnergy) {
    # interpolate the value between the two input arrays
    
    # +++ 20150804 I didn't think we should interpolate, not like described.  But what if tgtAzim is very close to azim_n, and far from azim_n-plus-1?
    #my $interpValue = ($inArray1[$iEnergy] + $inArray2[$iEnergy]) / 2;
    
    # +++ add comments
    # EA(theta_0, phi_1) = EA
    
    my $interpValue = $inArray1[$iEnergy] + 
                      ( ($varNaught-$var1) / ($var2-$var1) ) * 
                      ( $inArray2[$iEnergy] - $inArray1[$iEnergy]);
    
    # assign the interpolated value to the output
    $outArray[$iEnergy] = $interpValue;
  }

  ahlog::ah_debug "---- end interpolate ---- \n";
  
  # return a reference to the output array
  return (\@outArray);
  
} # end-sub interpolate

#############################################################################

# to compare floating point numbers
# equal(NUM1, NUM2, PRECISION) : returns true if NUM1 and NUM2 are
# equal to PRECISION number of decimal places
# from https://www.safaribooksonline.com/library/view/perl-cookbook/1565922433/ch02s03.html
sub equal {
    my ($A, $B, $dp) = @_;
    return sprintf("%.${dp}g", $A) eq sprintf("%.${dp}g", $B);
} # end-sub equal

#############################################################################

sub writeExtension {

  # -------------------------------------
  # define and initialize all variables
  # -------------------------------------
  
  # input parameters
  my $currFileName    = shift;
  my $currExtName     = shift;
  my $currOffaxis     = shift;
  my $currAzimuth     = shift;
  my $eaRef           = shift;
  my @ea              = @{$eaRef};
  
  # -------------------------------------
  
  ahlog::ah_debug "Creating output file $currFileName";
  
  ahlog::ah_debug "numUnqEnergies = $numUnqEnergies \n" ;
  
  # get the range of energy values.  currUnqEnergies is leftover from the last 
  # iteration through the history files, but it should have been the same for 
  # every file.
  my $energyValues = $currUnqEnergies[0] . "-" . $currUnqEnergies[$numUnqEnergies-1];

  # +++ 20150813 KLR change if vignetting
  # These are for CBD keywords
  my $cbdOffaxs   = "THETA(" . $currOffaxis . ")[arcmin]";
  my $cbdAzim     = "PHI(" . $currAzimuth . ")[deg]";
  my $cbdEnergy   = "ENERG(" . $energyValues . ")[keV]";
  
  # create text files to hold the column definitions, column data, and keywords
  # ftcreate needs these to create the output file
  my $columndeffile = "tmp_columns.txt";
  my $datafile      = "tmp_data.txt";
  my $headfile      = "tmp_head.txt";
  
  # add these to tmparray to cleanup later
  ahapp::add_temp_file($columndeffile);
  ahapp::add_temp_file($datafile);
  ahapp::add_temp_file($headfile);

  # the column name, unit, and comment will depend on if this is vignetting
  my $eaColumnName;
  my $eaColumnUnit;
  my $eaColumnComm;
  if ($doVignetting) {
    $eaColumnName = "EARATIO";
    $eaColumnUnit = "";
    $eaColumnComm = "Ratio of Effective Area per energy";
  } else {
    $eaColumnName = "EFFAREA";
    $eaColumnUnit = "cm**2";
    $eaColumnComm = "Effective Area per energy";
  }

  # create the column definition file, overwriting if it exists
  open my $fileHandleColumn, ">", $columndeffile or die "Can't open '$columndeffile'\n";
  print $fileHandleColumn "Energy 1D keV\n";
  print $fileHandleColumn "$eaColumnName 1D $eaColumnUnit\n";
#  if ($doVignetting) {
#    print $fileHandleColumn "EARATIO 1D\n";
#  } else {
#    print $fileHandleColumn "EFFAREA 1D cm**2\n";
#  }
  close $fileHandleColumn;
  
  # +++ 20150814 KLR for vignetting: change column name to EA_RATIO, no units
  
  # create the column data file, overwriting if it exists
  open my $fileHandleData, ">", $datafile or die "Can't open '$datafile'\n";
  for (my $iEnergy = 0 ; $iEnergy < $numUnqEnergies ; ++$iEnergy) {
    print $fileHandleData $currUnqEnergies[$iEnergy] . " " . $ea[$iEnergy] . "\n";
  }
  close $fileHandleData;

  # create the column definition file, overwriting if exists
  open my $fileHandleHeader, ">", $headfile or die "Can't open '$headfile'\n";
  print $fileHandleHeader "TELESCOP  $telescop        /    Mission name\n";
  print $fileHandleHeader "INSTRUME  $instrume        /    Instrument name\n";
  print $fileHandleHeader "DETNAM    $detnam          /    Telescope name\n";
  print $fileHandleHeader "CVSD0001  $cvsd0001        /    UTC date when file should be first used\n";
  print $fileHandleHeader "CVST0001  $cvst0001        /    UTC time when file should be first used\n";
  print $fileHandleHeader "CBD10001  $cbdOffaxs       /    Offaxis angle\n";
  print $fileHandleHeader "CBD20001  $cbdAzim         /    Azimuthal roll angle\n";
  print $fileHandleHeader "CBD30001  $cbdEnergy       /    Energy values\n";
  print $fileHandleHeader "OFFAXIS   $currOffaxis     /    [arcmin] Source off-axis angle\n";             # +++ 20150813 KLR change if vignetting
  print $fileHandleHeader "AZIMUTH   $currAzimuth     /    [deg] Source rotational (azimuthal) angle\n";  # +++ 20150813 KLR change if vignetting
  print $fileHandleHeader "Energy    $energyValues    /    [keV] Energy or energy range of image\n";
  print $fileHandleHeader "NOFFAXIS  $numUnqOffaxis   /    Number of off-axis angles\n";
  print $fileHandleHeader "NAZIMUTH  $numUnqAzimuth   /    Number of azimuthal angles\n";
  print $fileHandleHeader "NORMRAD   $normRadmm       /    [mm] Normalization radius\n";
  close $fileHandleHeader;

  # create the output file
  $status = ahgen::run_ftool("ftcreate",
                   "cdfile=$columndeffile",
                   "datafile=$datafile",
                   "outfile=$currFileName",
                   "headfile=$headfile",
                   "extname=$currExtName",
                   "clobber=yes");
  if ($status) {
    ahlog::ah_err "Error creating file $currFileName";
    ahapp::end_processing(1);
  }

  # add comments to the TTYPEn keywords
  $status = ahgen::run_ftool("fthedit",
                   "infile=$currFileName",
                   "keyword=TTYPE1",
                   "operation=add",
                   "value=Energy",
                   "comment=Energy");
  $status = ahgen::run_ftool("fthedit",
                   "infile=$currFileName",
                   "keyword=TTYPE2",
                   "operation=add",
                   "value=$eaColumnName",
                   "comment=$eaColumnComm");
  if ($status) {
    ahlog::ah_err "Error editing TTYPEn keywords in $currFileName";
    ahapp::end_processing(1);
  }

  # stamp parameters to output FITS as HISTORY keywords
  if ($ahapp::history) {   
    ahapp::write_parameters($currFileName, $currExtName) ;
  }
  
  # go to primary extension to write some additional keywords
  ahgen::set_keyword($currFileName, 0, "NOFFAXIS", 1 . " / Number of off-axis angles");
  ahgen::set_keyword($currFileName, 0, "NAZIMUTH", 1 . " / Number of azimuthal angles");
  ahgen::set_keyword($currFileName, 0, "NORMRAD", $normRadmm . " / [mm] Normalization radius");


  ahlog::ah_debug "---- end writeExtension ---- \n";
  
} # end sub writeExtension()

#############################################################################




# Revision Log:
# $Log: arftable.pl,v $
# Revision 1.9  2016/03/23 01:56:11  klrutkow
# add cvs log
#id:$ &
#

