#!/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) text file which lists the coordinates of regions of interest, one per line
#
# The region file should have the following format:
#
#       unit
#       x y inner_radius outer_radius 
#
#    where unit = unit the coordinates and radii are in; can be either d (decimal degrees) 
#                  or p (pixels) 
#          inner_radius = inner radius of annulus region to extract. Set this 
#                         to 0 to extract a circular region.
#          outer_radius = outer radius of annulus region to extract. This is 
#                         equal to the radius of a circular region. 
#

my $currDir;
chomp ($currDir = `pwd`);
my $unittype = 'N';
my $count = 0;
my $j;
my $currLine;
my @numbers;
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 $baseDir;
if ($currDir =~ /^(.*?[0-9]{10})/) {
    $baseDir = $1 ;
}

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 $position_file = <>);
$position_file = $currDir."/".$position_file;


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 @xval;
my @yval;
my @innerRad;
my @outerRad;
my @vals;
my $countlines=0;

#  Does the coord file have the proper format?

if (-e $position_file) { 
   print "Reading position file \n";
   open(COORDS, "<", $position_file); 
   while ($currLine = <COORDS>) {
      if ($currLine =~ m/D/i) {
         $unittype = 'D';              # decimal degrees
      }
      if ($currLine =~ m/P/i ) {
         $unittype = 'P';              # pixel
      }
      if ($unittype eq 'N') {
         print "No units found in the coords file. Exiting. \n";
         exit;
      }
      if ($currLine =~ m/\d/) {
         $countlines ++;
         # count numbers to confirm the correct format assign to variables
         @numbers = split (/\s/, $currLine);
         if (@numbers != 4) {
            print "The input file does not have the correct format. \n";
            exit;
         }
         push (@xval, $numbers[0]); 
         push (@yval, $numbers[1]); 
         push (@innerRad, $numbers[2]);
         push (@outerRad, $numbers[3]); 
      }
   }
   close (COORDS);
}
else {
   print "No region file found. \n" ;
   exit;
}

#
# If necessary, change the degree decimals to pixel values. We don't know
# which detectors the user wants, so we have to prepare for each case. 
#

my $convolved = 0;

if ($unittype eq 'D' ) { 
   if (($process_m1 == 1) and ($convolved == 0)) { 
      `evselect table=$m1_evt withimageset=yes imageset=m1_img_hold.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
      for ($j=0; $j<$countlines; $j++) {
         `ecoordconv imageset=m1_img_hold.fits x=$xval[$j] y=$yval[$j] coordtype=eqpos > hold.txt`; 
         open (HOLD, "<", "hold.txt");
         while ($currLine = <HOLD>)  {
            if ($currLine =~ m/ X: Y: / ) {
               @vals = split (/\s/, $currLine);
               $xval[$j] = $vals[3];
               $yval[$j] = $vals[4];
            }
         }
         close (HOLD);
         $innerRad[$j] = $innerRad[$j]/1.38889e-5;
         $outerRad[$j] = $outerRad[$j]/1.38889e-5;
      }
      `rm hold.txt`;
      `rm m1_img_hold.fits`;
      $convolved = 1;
   } 
   if (($process_m2 == 1) and ($convolved == 0)) { 
      `evselect table=$m2_evt withimageset=yes imageset=m2_img_hold.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
      for ($j=0; $j<$countlines; $j++) {
         `ecoordconv imageset=m2_img_hold.fits x=$xval[$j] y=$yval[$j] coordtype=eqpos > hold.txt`; 
         open (HOLD, "<", "hold.txt");
         while ($currLine = <HOLD>)  {
            if ($currLine =~ m/ X: Y: / ) {
               @vals = split (/\s/, $currLine);
               $xval[$j] = $vals[3];
               $yval[$j] = $vals[4];
            }
         }
         close (HOLD);
         $innerRad[$j] = $innerRad[$j]/1.38889e-5;
         $outerRad[$j] = $outerRad[$j]/1.38889e-5;
      }
      `rm hold.txt`;
      `rm m2_img_hold.fits`;
      $convolved = 1;
   }
   if (($process_pn == 1) and ($convolved == 0)) {
      `evselect table=$pn_evt withimageset=yes imageset=pn_img_hold.fits xcolumn=X ycolumn=Y imagebinning=binSize ximagebinsize=32 yimagebinsize=32`;
      for ($j=0; $j<$countlines; $j++) {
         `ecoordconv imageset=pn_img_hold.fits x=$xval[$j] y=$yval[$j] coordtype=eqpos > hold.txt`;
         open (HOLD, "<", "hold.txt");
         while ($currLine = <HOLD>)  {
            if ($currLine =~ m/ X: Y: / ) {
               @vals = split (/\s/, $currLine);
               $xval[$j] = $vals[3];
               $yval[$j] = $vals[4];
            }
         }
         close (HOLD);
         $innerRad[$j] = $innerRad[$j]/1.38889e-5;
         $outerRad[$j] = $outerRad[$j]/1.38889e-5;
      }
      `rm hold.txt`;
      `rm pn_img_hold.fits`;
   }
}

# 
# Make the spectrum. 
# 
# If inner radius = 0, use a circle. If not, use an annulus.
#

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) {
   for ($j=0; $j<$countlines; $j++) { 
      if ($innerRad[$j] == 0 ) {
         `evselect table='m1_evt_cp' energycolumn='PI' withfilteredset=yes filteredset='m1_filtered_$j.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE($xval[$j],$yval[$j],$outerRad[$j]))' withspectrumset=yes spectrumset='m1_src_$j.pi' spectralbinsize=15  withspecranges=yes specchannelmin=0 specchannelmax=11999 ` ; 
      } else {
         `evselect table='m1_evt_cp' energycolumn='PI' withfilteredset=yes filteredset='m1_filtered_$j.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE($xval[$j],$yval[$j],$outerRad[$j]))&&!((X,Y) IN circle($xval[$j],$yval[$j],$innerRad[$j]))' withspectrumset=yes spectrumset='m1_src_$j.pi' spectralbinsize=15  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`;
      `arfgen arfset=m1_src_$j.arf spectrumset=m1_src_$j.pi withrmfset=yes rmfset=m1_src_$j.rmf withbadpixcorr=yes badpixlocation=m1_evt_cp`; 
   }
}

if ($process_m2) {   
   for ($j=0; $j<$countlines; $j++) {
      if ($innerRad[$j] == 0 ) { 
         `evselect table='m2_evt_cp' energycolumn='PI' withfilteredset=yes filteredset='m2_filtered_$j.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE($xval[$j],$yval[$j],$outerRad[$j]))' withspectrumset=yes spectrumset='m2_src_$j.pi' spectralbinsize=15 withspecranges=yes specchannelmin=0 specchannelmax=11999 ` ;
      } else {
         `evselect table='m2_evt_cp' energycolumn='PI' withfilteredset=yes filteredset='m2_filtered_$j.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE($xval[$j],$yval[$j],$outerRad[$j]))&&!((X,Y) IN circle($xval[$j],$yval[$j],$innerRad[$j]))' withspectrumset=yes spectrumset='m2_src_$j.pi' spectralbinsize=15  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`;
      `arfgen arfset=m2_src_$j.arf spectrumset=m2_src_$j.pi withrmfset=yes rmfset=m2_src_$j.rmf withbadpixcorr=yes badpixlocation=m2_evt_cp`;
   }
}

if ($process_pn) {
   for ($j=0; $j<$countlines; $j++) {
      if ($innerRad[$j] == 0 ) {
         `evselect table='pn_evt_cp' energycolumn='PI' withfilteredset=yes filteredset='pn_filtered_$j.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE($xval[$j],$yval[$j],$outerRad[$j]))' withspectrumset=yes spectrumset='pn_src_$j.pi' spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=20479`; 
      } else { 
      `evselect table='pn_evt_cp' energycolumn='PI' withfilteredset=yes filteredset='pn_filtered_$j.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE($xval[$j],$yval[$j],$outerRad[$j]))&&!((X,Y) in CIRCLE($xval[$j],$yval[$j],$innerRad[$j]))' 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`;
      `arfgen arfset=pn_src_$j.arf spectrumset=pn_src_$j.pi withrmfset=yes rmfset=pn_src_$j.rmf withbadpixcorr=yes badpixlocation=pn_evt_cp`;
   } 
}

# 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";