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

# 
# Makes exposure maps for event files supplied by user, with 
# the low and high energies supplied in a text file. If no energy file exists, 
# the SAS defaults will be used.
#
# Requirements:
#    0) must be in same directory as input event files. 
#    1) text file which lists the events files to process,
#       one per line
#    2) If energy file exists, it should have the following format: 
#
#       lowE1 highE1
#       lowE2 highE2
#        ...   ...
#       lowEn highEn
#
#    where energies are in units of eV.
#

my $currDir;
chomp ($currDir = `pwd`);
my $energy_file="filler";
my $currLine;
my @numbers;
my @files;
my $file;
my $m1_evt;
my $m2_evt;
my $pn_evt;

my $process_pn = 0;
my $process_m1 = 0;
my $process_m2 = 0;
my $inst;
my $energyflag;
my $status = 0; 

my $odfDir;
my $test = 0;
 
# 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";
}

# On to the main event.

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

print "Should I use an energy file? If so, enter the file name; otherwise, enter N. ";
chomp (my $ans = <>) ;
if (($ans eq "N") or ($ans eq "n")) {
   $energyflag = 0 ;
} else {
   $energyflag = 1;
   $energy_file = $currDir."/".$ans;
}


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; stopping here. \n" ;
   exit;
}

#
# Does the energy file exist and have proper format?
#

my $countlines = 0;
my @loEn; 
my @hiEn;
my @loEnkeV;
my @hiEnkeV;
my $j=0;

if ($energyflag == 1) {
   if (-e $energy_file) {
      open(ENERGY, "<", $energy_file);
      while ($currLine = <ENERGY>) {
      $countlines ++;
      @numbers = split (/\s/, $currLine);
      # count array elements to confirm the correct format and assign to variables
      if (@numbers != 2) { 
         print "The input file does not have the correct format. \n";
         exit;
      }
      push (@loEn, $numbers[0]);
      push (@hiEn, $numbers[1]);
      push (@loEnkeV, $numbers[0]/1000.);
      push (@hiEnkeV, $numbers[1]/1000.);
      }
      close(ENERGY);
   }
   else {
      print "No energy file found; stopping here. \n";
      exit;
   }
}

#
# Trying to make this as non-fussy as possible, so even though the 
# user might have already produced images, we're going to make them 
# again anyway. Shouldn't add that much time to processing. 
#

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

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

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


chdir $currDir.'/expmap';

`atthkgen atthkset=atthk.fits timestep=1`;
if (! $energyflag) {
   if ($process_m1) {
      print "Making the M1 exposure map... \n";
      `evselect table=m1_evt_cp withimageset=yes imageset=m1_img_hold.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
      `eexpmap imageset=m1_img_hold.fits attitudeset=atthk.fits eventset=m1_evt_cp withdetcoords=false withvignetting=true usefastpixelization=true pimin=2000 pimax=4500 expimageset=m1_expmap_2-4.5keV.fits`;
      `rm m1_evt_cp m1_img_hold.fits`; 
   }
   if ($process_m2) {
      print "Making the M2 exposure map... \n";
      `evselect table=m2_evt_cp withimageset=yes imageset=m2_img_hold.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
      `eexpmap imageset=m2_img_hold.fits attitudeset=atthk.fits eventset=m2_evt_cp withdetcoords=false withvignetting=true usefastpixelization=true pimin=2000 pimax=4500 expimageset=m2_expmap_2-4.5keV.fits`;
      `rm m2_evt_cp m2_img_hold.fits`; 
   } 
   if ($process_pn) { 
      print "Making the PN exposure map... \n";
      `evselect table=pn_evt_cp withimageset=yes imageset=pn_img_hold.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
      `eexpmap imageset=pn_img_hold.fits attitudeset=atthk.fits eventset=pn_evt_cp withdetcoords=false withvignetting=true usefastpixelization=true pimin=2000 pimax=4500 expimageset=pn_expmap_2-4.5keV.fits`;
      `rm pn_evt_cp pn_img_hold.fits`;
   }
}

if ($energyflag eq 1) { 
   if ($process_m1 eq 1) {
      print "Making the M1 exposure map... \n";
      `evselect table=m1_evt_cp withimageset=yes imageset=m1_img_hold.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
      for ($j=0; $j<$countlines; $j++) {
      `eexpmap imageset=m1_img_hold.fits attitudeset=atthk.fits eventset=m1_evt_cp withdetcoords=false withvignetting=true usefastpixelization=true pimin=$loEn[$j] pimax=$hiEn[$j] expimageset=m1_expmap_$loEnkeV[$j]-$hiEnkeV[$j]keV.fits`;
      }
      `rm m1_evt_cp m1_img_hold.fits`;
   }
   if ($process_m2 eq 1) {
      print "Making the M2 exposure map... \n";
      `evselect table=m2_evt_cp withimageset=yes imageset=m2_img_hold.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
      for ($j=0; $j<$countlines; $j++) {
      `eexpmap imageset=m2_img_hold.fits attitudeset=atthk.fits eventset=m2_evt_cp withdetcoords=false withvignetting=true usefastpixelization=true pimin=$loEn[$j] pimax=$hiEn[$j] expimageset=m2_expmap_$loEnkeV[$j]-$hiEnkeV[$j]keV.fits`;
       }
      `rm m2_evt_cp m2_img_hold.fits`;
   }
   if ($process_pn eq 1) {
      print "Making the PN exposure map... This may take a while. Why don't you get up and get some coffee? \n";
      `evselect table=pn_evt_cp withimageset=yes imageset=pn_img_hold.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
      for ($j=0; $j<$countlines; $j++) {
      `eexpmap imageset=pn_img_hold.fits attitudeset=atthk.fits eventset=pn_evt_cp withdetcoords=false withvignetting=true usefastpixelization=true pimin=$loEn[$j] pimax=$hiEn[$j] expimageset=pn_expmap_$loEnkeV[$j]-$hiEnkeV[$j]keV.fits`;
       }
       `rm pn_evt_cp pn_img_hold.fits`;
   }
}

print "Done. \n";