#!/usr/bin/perl
use strict;
use warnings;

# 
# Makes light curve over entire FOV for event files supplied by user. 
# If a region file is also given, the light curve will be extracted 
# over both the full FOV and the region. 
#
#
# Requirements:
#    0) must be in same directory as input event files. 
#    1) text file which lists the events files to process, one per line
#    
# Optional: 
#    0) region file which lists the regions of interest, one per line, in phyical units. 
#       Must be in ds9 region file format. SAS can accept only circles, ellipses, boxes, 
#       and annuli.
#

my $currDir;
chomp ($currDir = `pwd`);
my $unittype = 'N';
my $count = 0;
my $currLine;
my @extraction_region;
my $source;
my $num;
my @files;
my $file;
my $m1_evt; 
my $m2_evt; 
my $pn_evt; 
my $j;
my $time_bin;
my $status;
my $inst;

my $process_pn = 0;
my $process_m1 = 0;
my $process_m2 = 0;
my $lcbin = 0 ; 
my $region_file;

my $region_flag;

####################################
# Set the environment variables.
####################################

my @res1;
my $baseDir;

if ($currDir =~ m/proc/) {
   @res1 = split(/proc/, $currDir);
   $baseDir = $res1[0];
}

if (-d $baseDir."/ODF") { 
   $ENV{SAS_ODFPATH} = $baseDir."/ODF";
   my $sumsas = glob ($baseDir."/ODF/*SUM.SAS");
   $ENV{SAS_ODF} = $sumsas;
   $ENV{SAS_CCF} = $baseDir."/ODF/ccf.cif";
} 
if (-d $baseDir."/odf") {
   $ENV{SAS_ODFPATH} = $baseDir."/odf";
   my $sumsas = glob ($baseDir."/odf/*SUM.SAS");
   $ENV{SAS_ODF} = $sumsas;
   $ENV{SAS_CCF} = $baseDir."/odf/ccf.cif";
}

####################################
# Get the event files.
####################################

print "Enter the file listing the event files: ";
chomp (my $event_file = <> );
$event_file = $currDir."/".$event_file;

print "Should I use a region file? If so, enter the file name; otherwise, enter N: ";
chomp (my $ans = <>) ;
if (($ans eq "N") or ($ans eq "n")) {
   $region_flag = 0 ;
} else { 
   $region_flag = 1;
   $region_file = $currDir."/".$ans;
   if ( ! -e $region_file) {
      print "There is no region file with that name. Stopping. \n" ;
      exit;
   }
}

print "Please enter the size of the time bins in seconds. ";
chomp ($lcbin = <>) ;
if ($lcbin !~ m/\d/) { 
   print "Time bin is not a number, stopping here. \n" ;
   exit;
}


if (-e $event_file) {
   open(EVENTS, "<", $event_file);
   while ($currLine = <EVENTS>) {
      $currLine = $currDir."/".$currLine;
      chomp ($currLine);
      @files = split (/\s/, $currLine);
      foreach $file (@files) {  # Check the headers to find out which instruments we need.
         $inst = $currDir."/"."instru.list";
         `fkeyprint $file instrume outfile=$inst`;
         open (INSTRU, "<", $inst);
         while ($currLine = <INSTRU>) {
            if ($currLine =~ m/INSTRUME= 'EPN     '/ ) {
               $pn_evt = $file;
               $process_pn = 1;
            }
            if ($currLine =~ m/INSTRUME= 'EMOS1   '/ ) {
               $m1_evt = $file;
               $process_m1 = 1;
            }
            if ($currLine =~ m/INSTRUME= 'EMOS2   '/ ) {
               $m2_evt = $file;
               $process_m2 = 1;
            }
         }
         close (INSTRU);
         `rm $inst`;
      }
   }
   close (EVENTS);
}
else {
   print "No events file found. \n" ;
   exit;
}

#######################################################
# Does the region file exist and have proper format?
# If so, grab the region info.
#######################################################

my $count_regions=0;
my $test_ds9 = 0;
my $test_physical = 0;

if ($region_flag) {
   if (-e $region_file) {
      print "Reading region file \n";
      open(REGION, "<", $region_file);
      while ($currLine = <REGION>) {
         if ($currLine =~ m/ DS9 /i) {
            $test_ds9 = 1;
         }
         if ($currLine =~ m'physical\n') {
            $test_physical = 1;
         }
      }
      close (REGION);
   }
   else {
      print "No region file found. \n" ;
      exit;
   }

   if ($test_ds9 == 0) {
      print "Region file is not in DS9 format. \n";
      exit;
   }
   if ($test_physical == 0) {
      print "Region file is not in correct units. \n";
      exit;
   }

   open(REGION, "<", $region_file);
   for ($j=0; $j<3; $j++) {   # Skip the header
      readline(REGION);
   }
   while ($currLine = <REGION>) {
      if (($currLine =~ m/circle/i) or ($currLine =~ m/ellipse/i) or ($currLine =~ m/box/i) or ($currLine =~ m/polygon/i) or ($currLine =~ m/annulus/i) or ($currLine =~ m/panda/i) or ($currLine =~ m/epanda/i) or ($currLine =~ m/bpanda/i)) {
         chomp($currLine);
         push (@extraction_region, $currLine);
         $count_regions ++;
      }
      else {
         print "Invalid shape in the region file! \n";
      }
   }
   close (REGION);
}

########################################
# Make the light curves of the full FOV.
########################################

if (! -d $currDir.'/lightcurve') { 
   `mkdir lightcurve`; 
} 

my $m1_evt_cp;
my $m2_evt_cp;
my $pn_evt_cp;

if ($process_m1 eq 1) { 
   $m1_evt_cp = $currDir."/lightcurve/m1_evt_cp";
   `cp $m1_evt $m1_evt_cp`;
}
if ($process_m2 eq 1) {
   $m2_evt_cp = $currDir."/lightcurve/m2_evt_cp";
   `cp $m2_evt $m2_evt_cp`;
}
if ($process_pn eq 1) {
   $pn_evt_cp = $currDir."/lightcurve/pn_evt_cp";
   `cp $pn_evt $pn_evt_cp`;
}

chdir $currDir.'/lightcurve';

if ($process_m1 eq 1) {
   `evselect table=m1_evt_cp withrateset=yes rateset='m1_ltcrv_fov_bin$lcbin.fits' maketimecolumn=yes timecolumn=TIME timebinsize=$lcbin makeratecolumn=yes`;
}
if ($process_m2 eq 1) {
   `evselect table=m2_evt_cp withrateset=yes rateset=m2_ltcrv_fov_bin$lcbin.fits maketimecolumn=yes timecolumn=TIME timebinsize=$lcbin makeratecolumn=yes`;
} 
if ($process_pn eq 1) { 
   `evselect table=pn_evt_cp withrateset=yes rateset=pn_ltcrv_fov_bin$lcbin.fits maketimecolumn=yes timecolumn=TIME timebinsize=$lcbin makeratecolumn=yes`;
}

##################################################################
# Make the light curves of the sources, if there is a region file.
##################################################################

if ($region_flag eq 1) {   
   for ($j=0; $j<$count_regions; $j++) {
      $source = shift(@extraction_region);
      print "Making light curves of the source(s)... \n";
      if ($process_m1 eq 1) {
         for ($j=0; $j<$count_regions; $j++)  {
            `evselect table=m1_evt_cp withrateset=yes rateset="m1_ltcrv_src$j\_bin$lcbin.fits" maketimecolumn=yes timecolumn=TIME timebinsize=$lcbin makeratecolumn=yes expression='(X,Y) IN $source'`;
         }
      }
      if ($process_m2 eq 1) {
         for ($j=0; $j<$count_regions; $j++)  {
            `evselect table=m2_evt_cp withrateset=yes rateset=m2_ltcrv_src$j\_bin$lcbin.fits maketimecolumn=yes timecolumn=TIME timebinsize=$lcbin makeratecolumn=yes expression='(X,Y) IN $source'`;
         }
      }
      if ($process_pn eq 1) {
         for ($j=0; $j<$count_regions; $j++)  {
            `evselect table=pn_evt_cp withrateset=yes rateset=pn_ltcrv_src$j\_bin$lcbin.fits maketimecolumn=yes timecolumn=TIME timebinsize=$lcbin makeratecolumn=yes expression='(X,Y) IN $source' `;
         }
      }
   }
}

# Let's be tidy now.

if ($process_m1) { 
   `rm m1_evt_cp`;
}
if ($process_m2) {
   `rm m2_evt_cp`;
}
if ($process_pn) {
   `rm pn_evt_cp`;
}

print "Done. \n";