#!/usr/local/bin/perl
#===============================================================================
#
# FILE NAME:    @(#)create_exposure_map	1.7
#
# DESCRIPTION:  A script to create an ACIS exposure map.  The default operation
#               of this script is to prompt the user before overwriting a file.
#               This behavior can be overridden by using the --clobber flag.
#
# AUTHOR:       Scott Koch (tsk@astro.psu.edu)
#               Copyright (C) 2000, Pennsylvania State University
# 
# NOTES:        
#
#               To view the aspect offsets data in IDL, use these commands:
#
#               a = mrdfits('aspect_offsets.fits',1,COLUMNS=['TIME', $
#                           'X_OFFSETS','Y_OFFSETS','ROLL_OFFSETS'])
#               function_1d,id,a.TIME,a.X_OFFSETS
#               function_1d,id,a.TIME,a.Y_OFFSETS
#               function_1d,id,a.TIME,a.ROLL_OFFSETS
#
#               Use these IDL commands to view the aspect histogram:
#
#               e = mrdfits('asphist_0.fits',1)
#               help,e,/st
#               dataset_2d,id,e.X_OFFSET,e.Y_OFFSET,weight=e.DURATION, $
#                          XBIN=1,YBIN=1
#
# REVISION HISTORY:
#
# July 24, 2000 - tsk - Handle EXPOSURE[0-9] keywords in a completely different
#                       way.  Instead of summing total CCD exposure then 
#                       weighting the exposure map by that, simply multiply
#                       exposure map for ccd N by the keyword EXPOSUR{N}.
#
# Aug 27, 2000  - tsk - Add a prompt so user can indicate whether to use a
#                       spectrum file with mkinstmap.  The default is NONE.
#===============================================================================
sub USAGE {
    printf  "\n\nUsage:  create_exposure_map {--clobber} xmin xmax ymin ymax binsize event_list aspect_file\n\n";
}
#===============================================================================
#
use File::Basename;
use Getopt::Long;

$result = GetOptions("clobber");

#-------------------------------------------------------------------------------
# Extract command line arguments
#-------------------------------------------------------------------------------
$xmin    = $ARGV[0] or die USAGE();
$xmax    = $ARGV[1] or die USAGE();
$ymin    = $ARGV[2] or die USAGE();
$ymax    = $ARGV[3] or die USAGE();
$binsize = $ARGV[4] or die USAGE();

#-------------------------------------------------------------------------------
# Check that the event file exists
#-------------------------------------------------------------------------------
$event_file = $ARGV[5] or die USAGE();
unless(-e $event_file) {die "\nFile not found:  $event_file\n\n";}


#-------------------------------------------------------------------------------
# Check that the aspect offset file exists
#-------------------------------------------------------------------------------
$aspect_file = $ARGV[6] or die USAGE();
unless(-e $aspect_file) {die "\nFile not found:  $aspect_file\n\n";}

#-------------------------------------------------------------------------------
# Use the pget command in the FTOOLS package, **NOT** the one in the ASCDS.
#-------------------------------------------------------------------------------
($name,$dir,$suffix) = fileparse(`which fdump`);
$pget = "${dir}pget";
unless (-e $pget) {die "Cannot find FTOOLS version of pget\n"};
$pset = "${dir}pget";
unless (-e $pset) {die "Cannot find FTOOLS version of pset\n"};

#-------------------------------------------------------------------------------
# determine which CCDs are in this dataset
#-------------------------------------------------------------------------------
SpawnCommand("fkeypar ${event_file}" . "+1 DETNAM");

# Parse the detector keyword and extract the active CCD_IDs for this
# observation.  The value should be something like 'ACIS-abcd' (with single
# quotes).
@result = SpawnCommand("$pget fkeypar value");
$detnam = $result[0];
chomp($detnam);
$active_ccds = [];
for ($ii=6; $ii<length($detnam)-1; $ii++) {
  $ccd = substr($detnam,$ii,1);
  unless ( ($ccd >= 0) && ($ccd <=9) )  {
    die "\nError extracting active CCDs from DETNAM keyword:  $detnam\n\n";
  }
  @active_ccds = (@active_ccds, $ccd);
}
$num_ccds = (scalar @active_ccds);

printf "Active CCDs are:";
for ($ii=0; $ii<$num_ccds; $ii++) {
  printf "  $active_ccds[$ii]";
}
printf "\n";

#-------------------------------------------------------------------------------
# Fetch a filter string to get rid of periods with bad aspect offsets.
#-------------------------------------------------------------------------------
printf "\nEnter a criterion for filtering the aspect offsets file:  \n";
$aspect_filter = <STDIN>;
chomp ($aspect_filter);

#-------------------------------------------------------------------------------
# Prompt user for use of a spectrum file.
#-------------------------------------------------------------------------------
printf "\nEnter the name of an appropriate spectrum file (used by mkinstmap).  Press RETURN to use the default energy:  \n";
$spectrumfile = <STDIN>;
chomp ($spectrumfile);
if($spectrumfile eq "") {
  $spectrumfile = "NONE";
} else {
  unless(-e $spectrumfile) {die "\nFile not found:  $spectrumfile\n\n";}
}

#-------------------------------------------------------------------------------
# Create a gti file which only includes periods of "good" aspect solutions.
#-------------------------------------------------------------------------------
printf "\nFiltering aspect offset file by above criterion...\n\n";
$gtifile = 'gti.fits';
if(OverwriteFile($gtifile)) {
  SpawnCommand("dmgti infile=$aspect_file outfile=$gtifile userlimit=\'$aspect_filter\' clobber='yes'");
}
# Due to a bug in dmgti, we have to subtract a very small amount of time off of
# all of the STOP times.
$fixed_gtifile = 'gti.fixed.fits';
if(OverwriteFile($fixed_gtifile)) {
  SpawnCommand("fcalc infile=$gtifile+2 outfile=$fixed_gtifile clname=\"STOP\" expr=\"STOP-0.00001\" clobber=yes");
}

#-------------------------------------------------------------------------------
# Filter the original aspect offsets file with the "good" gti file 
#-------------------------------------------------------------------------------
printf "\nFiltering the Aspect Offsets file, $aspect_file, by the gti file, $gtifile ...\n\n";
$filtered_aspect_file = 'aspect_offsets.fits';
if(OverwriteFile($filtered_aspect_file)) {
  SpawnCommand("dmcopy infile=\'$aspect_file\[\@$fixed_gtifile\]\' outfile=$filtered_aspect_file clobber=yes");
}

#-------------------------------------------------------------------------------
# Now filter the input event list by the gti table.  This should get rid of
# the data for times when the aspect offsets are bad.
#-------------------------------------------------------------------------------
printf "\nFiltering input event list by gti file to remove data with bad aspect...\n\n";
$filtered_event_file = 'filtered_event_list.fits';
if(OverwriteFile($filtered_event_file)) {
  SpawnCommand("dmcopy infile=\'$event_file\[\@$fixed_gtifile\]\' outfile=$filtered_event_file clobber=yes");
}

@result = `fkeyprint ${event_file}[EVENTS] NAXIS2 | grep =`;
print "  Original event count in $event_file:  ",  @result, "\n";

@result = `fkeyprint ${filtered_event_file}[EVENTS] NAXIS2 | grep =`;
print "\n  Filtered event count in $filtered_event_file:  ", @result, "\n";
print "\n";

#-------------------------------------------------------------------------------
# Now create an image of the sky binned the same as the exposure map
#-------------------------------------------------------------------------------
$sky_image = "sky_image.fits";
printf "\nCreating a sky image...\n";
if(OverwriteFile($sky_image)) {
  SpawnCommand("dmcopy \"$filtered_event_file\[bin x=$xmin:$xmax:$binsize,y=$ymin:$ymax:$binsize\]\" $sky_image clobber=yes");
}

#-------------------------------------------------------------------------------
# Determine the number of pixels in the X and Y directions.
#------------------------------------------------------------------------------
`fkeypar ${sky_image}\[PRIMARY\] NAXIS1`;
$num_x_pixels = `$pget fkeypar value`;
chomp $num_x_pixels;
`fkeypar ${sky_image}\[PRIMARY\] NAXIS2`;
$num_y_pixels = `$pget fkeypar value`;
chomp $num_y_pixels;

#-------------------------------------------------------------------------------
# Apply SIM corrections.
#-------------------------------------------------------------------------------
$simfile = 'sim.fits';
if(OverwriteFile($simfile)) {
  printf "\nApplying SIM corrections to filtered Aspect Offsets file...\n\n";
  SpawnCommand("asp_apply_sim infile=$filtered_aspect_file outfile=$simfile clobber=yes");
}

#-------------------------------------------------------------------------------
# Create an aspect histogram file. Make an instrument map and exposure map.
#-------------------------------------------------------------------------------
@exposure_times = ();
for ($ii=0; $ii<$num_ccds; $ii++) {

  printf "================================================================================\n";
  printf "Creating exposure map for CCD %d\n", $active_ccds[$ii];

  $ccd = $active_ccds[$ii];

#-------------------------------------------------------------------------------
  printf "  Creating aspect histogram file ...\n";
#-------------------------------------------------------------------------------
  $outfile = "asphist_$ccd.fits";
  if(OverwriteFile($outfile)) {
    SpawnCommand("asphist infile=$simfile gtifile=\'$filtered_event_file\[ccd_id=$ccd\]\' outfile=$outfile dtffile=\'\' clobber=yes");
  }

#-------------------------------------------------------------------------------
  printf "  Creating instrument map file ...\n";
#-------------------------------------------------------------------------------
  $outfile = "instmap_$ccd.fits";
  if(OverwriteFile($outfile)) {
    SpawnCommand("mkinstmap obsfile=\'asphist_${ccd}.fits\[asphist\]\' outfile=$outfile det=ACIS-${ccd} mirror=HRMA minxpixel=1 minypixel=1 numxpixels=-1 numypixels=-1 numinstmapxpixels=1024 numinstmapypixels=1024 ardlibparfile=ardlib.par spectrumfile=${spectrumfile} monoenergy=4 verbose=0");
  }

#-------------------------------------------------------------------------------
  printf "  Creating exposure map file ...\n";
#-------------------------------------------------------------------------------
  $outfile = "expmap_${ccd}.fits";
  if(OverwriteFile($outfile)) {
    $num_sky_x = $xmax - $xmin;
    $num_sky_y = $ymax - $ymin;
    SpawnCommand("mkexpmap instmapfile=instmap_${ccd}.fits outfile=$outfile minskyx=$xmin minskyy=$ymin numskyx=$num_sky_x numskyy=$num_sky_y numexpmapxpixels=$num_x_pixels numexpmapypixels=$num_y_pixels asphistfile=\'asphist_${ccd}.fits\[asphist\]\' useavgaspect=no verbose=0");
  }

#-------------------------------------------------------------------------------
# Grab the exposure time from the event file.
#-------------------------------------------------------------------------------
  `fkeypar $filtered_event_file\[1\] EXPOSUR${ccd}`;
  $exposure_times[$ii] = `$pget fkeypar value`;
  printf "  Exposure time is %f\n", $exposure_times[$ii];

#-------------------------------------------------------------------------------
# Multiply the exposure map (in units of cm^2) by the exposure to yield an
# exposure map which has been corrected for effective area.  The resulting
# image is in units of (s cm^2).
#-------------------------------------------------------------------------------
  printf "Creating weighted exposure map for CCD %d\n", $active_ccds[$ii];
#-------------------------------------------------------------------------------
  $outfile = "wexmap_${ccd}.fits";
  if(OverwriteFile($outfile)) {
    SpawnCommand("dmimgcalc infile=expmap_${ccd}.fits infile2=expmap_${ccd}.fits outfile=$outfile operation=add weight=$exposure_times[$ii] weight2=0.0 verbose=1 clobber=yes");
  }
}

#-------------------------------------------------------------------------------
# Now combine the weighted exposure maps into a single binned exposure map
# image.
#-------------------------------------------------------------------------------
$expmap_all = "expmap_all.fits";
printf "\nCombining individual exposure maps into a single exposure map...\n";
`ls wexmap_[0-9].fits > wexmap.lis`;
if(OverwriteFile($expmap_all)) {
 SpawnCommand("dmregrid \@wexmap.lis $expmap_all rotang=0 npts=1 bin='1:$num_x_pixels:1,1:$num_y_pixels:1' xoffset=0 yoffset=0 rotxcenter=0 rotycenter=0 clobber=yes");
}

#-------------------------------------------------------------------------------
# Finally, divide the sky image by the exposure map.
#-------------------------------------------------------------------------------
$norm_image = "sky_image_norm.fits";
printf "\nNormalizing $sky_image by $expmap_all...\n";
if(OverwriteFile($norm_image)) {
  SpawnCommand("dmimgcalc infile=$sky_image infile2=$expmap_all outfile=$norm_image operation=div weight=1.0 weight2=1.0 verbose=1 clobber=yes");
}

printf "\ncreate_exposure_map complete.\n\n";
#===============================================================================
sub OverwriteFile {

  local $ret    = 1;
  local $file   = $_[0];
  local $answer = "";

  if ($opt_clobber) {
    return 1;
  }

  if (-e $file) {
    printf "$file exists.  Press 1 (one) to overwrite:  ";
    $answer = <STDIN>;
    chomp $answer;

    if($answer != 1) {
      $ret = 0;
    }
  }

  return $ret;
}

#===============================================================================
sub SpawnCommand {
  local $command = $_[0];

  printf "  $command\n";
  return `$command`;
}