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

#
# Makes the spectrum and response files for a given set of coordinates. 
# Requirements: 
#    0) ccf.cif and *SUM.SAS files exist and environment variables are set
#    1) run from proc directory
#    2) text file which lists the events files to process, one per line
#    3) region file which lists the regions of interest, one per line, in phyical units
#
# The region file should be in the format of a ds9 region file. SAS can accept
# only circles, ellipses, boxes, and annuli. 
#

my $currDir;
chomp ($currDir = `pwd`);
my $unittype = 'N';
my $count = 0;
my $j;
my $currLine;
my $extended_flag;
my @extraction_region;
my @extraction_region_m1;
my @extraction_region_m2;
my @extraction_region_pn;
my $source;
my $num;
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;

# 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 "Enter the region file name: ";
chomp (my $region_file = <>);
$region_file = $currDir."/".$region_file;

print "Is this an extended source? y/n: ";
chomp (my $ans = <> );
if (($ans eq "N") or ($ans eq "n")) {
   $extended_flag = 0 ;
}
if (($ans eq "Y") or ($ans eq "y")) {
   $extended_flag = 1 ;
}
if (($ans ne "Y") and ($ans ne "y") and ($ans ne "N") and ($ans ne "n")) {
   print "I don't know what you mean. Stopping. \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;
}

my $count_regions=0;

#  Does the region file have the proper format?

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

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 spectrum. 
# 

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

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

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

chdir $currDir.'/spectrum';

print "Making the spectrum... \n";

if ($process_m1) {
   @extraction_region_m1 = @extraction_region;
   for ($j=0; $j<$count_regions; $j++) { 
      $source = shift(@extraction_region_m1);
      `evselect table='m1_evt_cp' energycolumn='PI' withfilteredset=yes filteredset='m1_filtered_$j.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in $source)' withspectrumset=yes spectrumset='m1_src_$j.pi' spectralbinsize=5  withspecranges=yes specchannelmin=0 specchannelmax=11999 ` ;
      print "Testing M1 for pileup. Output file is m1_epat_$j.ps. \n";
      `epatplot set=m1_filtered_$j.fits plotfile=m1_epat_$j.ps useplotfile=yes`;
      print "Making the M1 response and ancillary files. This may take a while. \n"; 
      `rmfgen rmfset=m1_src_$j.rmf spectrumset=m1_src_$j.pi`;
      if ($extended_flag == 0) {
      `arfgen arfset=m1_src_$j.arf spectrumset=m1_src_$j.pi withrmfset=yes rmfset=m1_src_$j.rmf withbadpixcorr=yes badpixlocation=m1_evt_cp setbackscale=yes` ; 
      }
      if ($extended_flag == 1) {
      `arfgen arfset=m1_src_$j.arf spectrumset=m1_src_$j.pi withrmfset=yes rmfset=m1_src_$j.rmf withbadpixcorr=yes badpixlocation=m1_evt_cp extendedsource=yes detmaptype=flat setbackscale=yes` ;
      }
   }
}

if ($process_m2) {   
   @extraction_region_m2 = @extraction_region;
   for ($j=0; $j<$count_regions; $j++) {
      $source = shift(@extraction_region_m2);
      `evselect table='m2_evt_cp' energycolumn='PI' withfilteredset=yes filteredset='m2_filtered_$j.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in $source)' withspectrumset=yes spectrumset='m2_src_$j.pi' spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=11999 ` ;
      print "Testing M2 for pileup. Output file is m2_epat_$j.ps. \n";
      `epatplot set=m2_filtered_$j.fits plotfile=m2_epat_$j.ps useplotfile=yes`;
      print "Making the M2 response and ancillary files. This may take a while. \n";
      `rmfgen rmfset=m2_src_$)j.rmf spectrumset=m2_src_$j.pi`;
      if ($extended_flag == 0) {
      `arfgen arfset=m2_src_$j.arf spectrumset=m2_src_$j.pi withrmfset=yes rmfset=m2_src_$j.rmf withbadpixcorr=yes badpixlocation=m2_evt_cp setbackscale=yes`;
      }
      if ($extended_flag == 1) {
      `arfgen arfset=m2_src_$j.arf spectrumset=m2_src_$j.pi withrmfset=yes rmfset=m2_src_$j.rmf withbadpixcorr=yes badpixlocation=m2_evt_cp extendedsource=yes detmaptype=flat setbackscale=yes` ;
      }
   }
}

if ($process_pn) {
   @extraction_region_pn = @extraction_region;
   for ($j=0; $j<$count_regions; $j++) {
      $source = shift(@extraction_region_pn);
      `evselect table='pn_evt_cp' energycolumn='PI' withfilteredset=yes filteredset='pn_filtered_$j.fits' keepfilteroutput=yes filtertype='expression' expression='(PATTERN <= 4) && (FLAG == 0) && ((X,Y) in $source)' withspectrumset=yes spectrumset='pn_src_$j.pi' spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=20479`; 
      print "Testing PN for pileup. Output file is pn_epat_$j.ps. \n";
      `epatplot set=pn_filtered_$j.fits plotfile=pn_epat_$j.ps useplotfile=yes`;
      print "Making the PN response and ancillary files. This may take a while. \n";
      `rmfgen rmfset=pn_src_$j.rmf spectrumset=pn_src_$j.pi`;
      if ($extended_flag == 0) {
      `arfgen arfset=pn_src_$j.arf spectrumset=pn_src_$j.pi withrmfset=yes rmfset=pn_src_$j.rmf withbadpixcorr=yes badpixlocation=pn_evt_cp setbackscale=yes`;
      }
      if ($extended_flag == 1) {
      `arfgen arfset=pn_src_$j.arf spectrumset=pn_src_$j.pi withrmfset=yes rmfset=pn_src_$j.rmf withbadpixcorr=yes badpixlocation=pn_evt_cp extendedsource=yes detmaptype=flat setbackscale=yes` ;
      }
   } 
}

# Next to godliness.

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

print "Done. \n";