#!/usr/local/bin/perl
#===============================================================================
#
# FILE NAME:    @(#)extract_sources.plx	1.6
#
# DESCRIPTION:  A Perl script that creates a series of event lists corresponding
#		to a series of circular extraction regions.
#
# The extraction regions are specified in a simple ASCII table in the file 
# "catalog" with one of these formats:
#	x	y	radius
# OR
#	x	y	radius	source_name
#
# The fields x,y,&radius are in units of sky pixels.
#
# If the field "source_name" is omitted then the event list generated for that
# source will have a name derived from the source number.
# The source name is recorded in the FITS keyword OBJECT in the output event 
# lists.
#
# The option --suffix can be used to add a suffix (e.g. the obsid) to the
# source names. 
#
# Normally each extracted event list will be put in its own directory, but if
# the option --flat is used they will all go in the current directory.
#
# In addition, a ds9-compatible region file depicting the circular regions is
# created with the name "source_regions.ds9".
#
# If -background specified, the script also creates an event list named 
# "background.evt" containing the events that do NOT fall in ANY of the 
# circular regions.
#
# *** WARNING!!!  
#     If any of the output files already exist they are OVERWRITTEN!!! 
# *** WARNING!!!  

# AUTHOR:       Patrick Broos (patb@astro.psu.edu)
#               Copyright (C) 2000, Pennsylvania State University
#
#===============================================================================

sub USAGE {
    printf  "\n\nUsage:  extract_sources.plx event_file catalog -suffix=source_name_suffix -flat -background\n      -rmf rmf_file -asp0 asphist_file0 -asp1 asphist_file1 -asp2 asphist_file2 -asp3 asphist_file3 -asp4 asphist_file4 -asp5 asphist_file5 -asp6 asphist_file6 -asp7 asphist_file7 -asp8 asphist_file8 -asp9 asphist_file9\n\n";
}
#===============================================================================

use File::Copy;
use Getopt::Long;

#-------------------------------------------------------------------------------
# Extract command line arguments
#-------------------------------------------------------------------------------
GetOptions("suffix=s", \$suffix, "flat", "background", "rmf=s", "asp0=s", "asp1=s", "asp2=s", "asp3=s", "asp4=s", "asp5=s", "asp6=s", "asp7=s", "asp8=s", "asp9=s");

$event_file  = $ARGV[0] or die USAGE();
$catalog = $ARGV[1] or die USAGE();

#-------------------------------------------------------------------------------
# Prepare output files.
$region_file     = "source_regions.ds9";
$background_file = "background.evt";
$temp_file	 = "extract_sources.temp";

#if (-e $region_file)     {unlink ${region_file};}
#if ($opt_background && -e $background_file) {unlink ${background_file};}
if (-e $temp_file)       {unlink ${temp_file};}

if (-e $region_file)     {die "REMOVE ${region_file}\n";}
if ($opt_background && -e $background_file) {die "REMOVE ${background_file}\n";}
if (-e $temp_file)       {die "REMOVE ${temp_file}\n";}

open(SOURCE_FILE,  "${catalog}") || die "\ncannot open $catalog\n";
open(REGION_FILE, ">${region_file}") || die "\ncannot open $region_file\n";
print REGION_FILE "physical\n";

if ($opt_background)
  {copy( ${event_file}, ${temp_file} );}

#-------------------------------------------------------------------------------
# Process source file.
# Define field delimiters as spaces and tabs.
# We do NOT use \s because we do NOT want to include \n as a delimiter!
# We use * multiplier (zero or more) so that initial whitespace is optional.
$ws = "[ \t]*";

@ccd_list = (0,1,2,3,4,5,6,7,8,9);
$source_num = 0;
while (<SOURCE_FILE>) 
  {
  # Parse a line of the ASCII table.
  # We use S* instead of S+ to keep the matching "greedy", i.e. to NOT split
  # a single number into multiple tokens.
  ($xx, $yy, $radius, $source_name) = /$ws(\S*)$ws(\S*)$ws(\S*)$ws(\S*)/;
  
  # The parsing statement above could be replaced by:
  # s/^$ws//;     remove leading whitespace
  # ($xx, $yy, $radius, $source_name) = split(/$ws/); 
  
  # Create names for this source and for its event list.
  if ($source_name EQ "") 
    {$source_name = sprintf("%4.4d",$source_num);}
    
  printf( "\n\n%s\nSOURCE: %s at (%f, %f)\n", "-" x 80, $source_name, $xx, $yy );
    
  if ($opt_flat) 
    {$path = "";}
  else
    {
    $path = "$source_name/";
    if (! -e $path) 
      {mkdir( $path, 0777 ) || die "Cannot make directory: $!";}
    }

  # Build a region specifier and write to REGION_FILE.
  $region = sprintf("circle(%f,%f,%f)", $xx, $yy, $radius);
  print REGION_FILE $region, "\n";
  
  # Execute dmcopy command to extract region.
  $source_file = sprintf("%s%s%s.evt", $path, $source_name, $suffix);
  
  SpawnCommand("dmcopy \"${event_file}\[ sky=${region} \]\" ${source_file} clobber=yes");
  
  # Add source name to FITS header.
  SpawnCommand("dmhedit ${source_file} operation=add key=OBJECT value=${source_name}");
  
  # Execute dmcopy command to exclude region from background_file.
  if ($opt_background)
    {
    rename ${background_file}, ${temp_file};
    SpawnCommand("dmcopy \"${temp_file}\[ sky=field()-${region} \]\" ${background_file}");
    }


  # If desired, build 5 ARF's which sample the extraction circle.
  if ($opt_rmf)
    {
    @arfs_to_sum = ();
    foreach $id (@ccd_list)
      {
      # Run dmlist on a CCD-filtered event list.
      $_ = `dmlist "${source_file}\[ccd_id=$id\]" blocks`;
      
      # Parse out the number of rows in the EVENTS HDU.
      /[.\n]*EVENTS.+x *(\d+)/;
      $num_events = $1;
      
      if ($num_events GT 0)
        {
        if ($id LE 3)
          {$detsubsys = "ACIS-I" . $id;}
        else
          {$detsubsys = "ACIS-S" . ($id - 4);}
        
        $asphist_file = "opt_asp${id}";
        if (! -r ${$asphist_file}) 
          {die "\n ERROR: You must supply valid -asp${id} parameter.\n";}
        
        @xoffset = (0, 0, $radius, 0, -$radius);
        @yoffset = (0, $radius, 0, -$radius, 0);
        
        printf( "\n%d events on %s\n", $num_events, $detsubsys );
        print "GENERATING 5 ARFs to sample the extraction region ...\n";

        @arfs_to_average = ();
        for ($ii=0; $ii<scalar(@xoffset); $ii++)
          {
          $sourcepixelx = $xx + $xoffset[$ii];
          $sourcepixely = $yy + $yoffset[$ii];
          
          $arf_file = sprintf("%s%s%s.%d_%d.arf", $path, $source_name, $suffix, $id, $ii);
          push( @arfs_to_average, $arf_file );
        
          SpawnCommand("mkarf outfile=$arf_file rmffile=$opt_rmf asphistfile=\"${$asphist_file}\[asphist\]\" sourcepixelx=$sourcepixelx sourcepixely=$sourcepixely detsubsys=$detsubsys mirror=HRMA grating=NONE ardlibparfile=ardlib.par obsfile=\"${$asphist_file}\[asphist\]\" verbose=0");
          
          } # loop over ARF positions
          
        # Average the 5 ARFs we used to cover the extraction circle.
        $weight = 1.0 / scalar(@arfs_to_average);
        open( TEMP_FILE, ">$temp_file") || die "\ncannot open $temp_file\n";
        foreach (@arfs_to_average) {print TEMP_FILE "$_ $weight\n";}
        close(TEMP_FILE);
        
        $arf_file = sprintf("%s%s%s.%d.arf", $path, $source_name, $suffix, $id);
        push( @arfs_to_sum, $arf_file );
        
        print "\nAVERAGING ARFs that sample the extraction region ...\n";
        SpawnCommand("addarf list=\@${temp_file} out_ARF=\!${arf_file} chatter=0");
        
        } # ($num_events GT 0)
      } # foreach CCD

    
    # ADD the ARF's from multiple CCDs if necessary.
    $arf_file = sprintf("%s%s%s.arf", $path, $source_name, $suffix);
    
    if (scalar(@arfs_to_sum) EQ 1)
      {rename $arfs_to_sum[0], $arf_file;}
    else
      {
      print "\nADDING ARFs from multiple CCDs ...\n";
      open( TEMP_FILE, ">$temp_file") || die "\ncannot open $temp_file\n";
      foreach (@arfs_to_sum) {print TEMP_FILE "$_ 1.0\n";}
      close(TEMP_FILE);
    
      SpawnCommand("addarf list=\@${temp_file} out_ARF=\!${arf_file} chatter=0");
      }
    } # opt_rmf
 
  $source_num++;
  }

close(REGION_FILE);
close(SOURCE_FILE);
if (-e $temp_file)       {unlink ${temp_file};}
exit;

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

  print "  $command\n";
  print `$command`;
}