Subs::XrtGrbLc (version $)


package Subs::XrtGrbLc;
##############################################################################
#
# DESCRIPTION: Produce XRT GRB lightcurve, gif, and related files.  Runs
# xrtgrblc for PC and WT mode events; xrtproducts for PD/LR and PD/PU mode
# events.
#
# HISTORY:
# HISTORY: $Log: XrtGrbLc.pm,v $
# HISTORY: Revision 1.21  2016/11/04 12:56:36  apsop
# HISTORY: Do not remove UTCFINIT from attitude files since it is used by prefilter to accomodate Swift clock drift.
# HISTORY:
# HISTORY: Revision 1.20  2016/09/29 15:08:24  apsop
# HISTORY: Fix calculation of angular separation
# HISTORY:
# HISTORY: Revision 1.19  2015/09/28 15:28:55  apsop
# HISTORY: Updated to HEASoft 6.17 and applied XRT, UVOT, and clock CALDB patches. Modified XRT event file processing to avoid further processing after an error occurs.
# HISTORY:
# HISTORY: Revision 1.18  2015/05/22 13:59:33  apsop
# HISTORY: Updates for 3.16.08:
# HISTORY: * modify how xrtgrblc and xrtexpomap are run as requested by ASDC.
# HISTORY: * avoid E2 in UVOTReport.pm that occurs when trying to run mgtime without GTI.
# HISTORY: * update calibration clock file to v107.
# HISTORY:
# HISTORY: Revision 1.17  2014/08/15 09:33:47  apsop
# HISTORY: make_xrt_grb_lc: More informative email messages to SDC staff if can't
# HISTORY: findGRB in catalogs.  Log a different message if it was too soon to
# HISTORY: send the email.
# HISTORY:
# HISTORY: Revision 1.16  2014/04/21 21:09:36  apsop
# HISTORY: Generate light curve, pha, and arf files for XRT PD/LR and
# HISTORY: PD/PU modes.  New subs make_pd_mode_lc and concatLogs.
# HISTORY: Required some reorganization to allow processing PD modes
# HISTORY: potentially without PC or WT modes.  write_standard_keyword
# HISTORY: now logs the file it's doing.  Several minor and commentary
# HISTORY: improvements.
# HISTORY:
# HISTORY: Revision 1.15  2014/02/28 11:58:29  apsop
# HISTORY: Don't use the uat attitude file, at request of XRT Team.
# HISTORY: No longer sets ATTFLAG keywords (leaves that for SW0WrapUp).
# HISTORY: Check for final run using jobpar object rather than job_title parameter.
# HISTORY: Include trigger ID in subject of the "GRB not found in catalogs" email.
# HISTORY:
# HISTORY: Revision 1.14  2013/08/10 06:18:56  apsop
# HISTORY: Use FITSfile objects instead of system() calls and backticks to
# HISTORY: set/get the ANCRFILE and RESPFILE keywords (others already did).
# HISTORY: system() and backticks use the login shell environment, which is
# HISTORY: set up for an old headas and is incompatible with heasoft-6.13
# HISTORY: (eg. PFILES); FITSfile runs tools in the correct environment as set
# HISTORY: in sw0.par. Also added error handling on getting the RMF files, and
# HISTORY: list the xspec script.  Added comments, esp. subroutine headers.
# HISTORY:
# HISTORY: Revision 1.13  2013/04/03 22:01:26  apsop
# HISTORY: sub sendEmail: strip apostrophes from subject, so they don't break
# HISTORY: the constructed command line.
# HISTORY:
# HISTORY: Revision 1.12  2012/08/30 07:36:30  apsop
# HISTORY: Many more log entries to tell which GRB coordinates are used and
# HISTORY: which catalog they came from, as in UvotProduct.pm and XrtProduct.pm.
# HISTORY: Reworded the email for Optical/XRT coordinates mismatch so it's clearer
# HISTORY: that is a warning and what should be checked.  Note source of trigtime
# HISTORY: (although info is not currently used).  Lots of minor formatting changes.
# HISTORY:
# HISTORY: Revision 1.11  2012/02/15 01:38:29  apsop
# HISTORY: The /* following setplot splashpage in the xspec command list should
# HISTORY: have been on its own line.
# HISTORY:
# HISTORY: Revision 1.10  2012/02/07 22:44:27  apsop
# HISTORY: Added object and sequence to "OPTICAL GRB coordinates
# HISTORY: outside XRT coordinates" alert message.
# HISTORY:
# HISTORY: Revision 1.9  2012/02/04 08:05:52  apsop
# HISTORY: Added "/*" to setplot splashpage off command, supposedly this avoids
# HISTORY: certain problems if used with an older version of Xspec.
# HISTORY:
# HISTORY: Revision 1.8  2012/02/02 06:22:39  apsop
# HISTORY: Added "setplot splashpage off" to xspec script.
# HISTORY:
# HISTORY: Revision 1.7  2012/01/12 06:52:03  apsop
# HISTORY: Changes going to proc3.15.03
# HISTORY:
# HISTORY: Revision 1.5  2011/01/20 19:00:12  apsop
# HISTORY: Added code to:
# HISTORY: 	1 - use UVOT attitude file if available
# HISTORY: 	2 - obtain correct TRIGTIME
# HISTORY: 	3 - use correct burst RA and DEC
# HISTORY:
# HISTORY:
#
# VERSION: $Revision: 1.21 $
#
##############################################################################

use Subs::Sub;
use Util::GTIfile;
use Util::SWCatalogue;

use File::Path;
use File::Basename;
use Cwd;

@ISA = ("Subs::Sub");
use strict;

use constant PI => atan2(0, -1);


sub new {
    my $proto=shift;
    my $self=$proto->SUPER::new();

    $self->{DESCRIPTION}="Make XRT GRB Light Curve";

    return $self;
}

##################
# METHODS:
##################

sub body {
    my $self=shift;

    my $log     =$self->log();
    my $filename=$self->filename();
    my $jobpar	=$self->jobpar();
    my $procpar  = $self->procpar();



    #########################
    # make xrt grb light curve
    #########################
    $self->make_xrt_grb_lc();

} # end of body method


#############################################################################
#
# make_xrt_grb_lc: Produce xrt grb light curve files
#
# Parameters: none
#
#############################################################################

sub make_xrt_grb_lc {
  my $self=shift;

  my $log      = $self->log();
  my $filename = $self->filename();
  my $jobpar   = $self->jobpar();
  my $procpar  = $self->procpar();


  my $watchers = $procpar->read("watchers");
  my $object = $jobpar->read('object');


  my $seq = $jobpar->read('sequence');
  my $segnum = substr($seq, 8);
  $segnum = $segnum*1;
  if ($segnum > 0) {
    $log->entry("Segment number is greater than zero: $segnum, no light curves or spectral analysis made");
    return;
  }


  ###################################
  # CHECK IF WE HAVE ALL FILES NEEDED
  ###################################

  # Event files for each mode.  The PC and WT modes can be processed
  # together, but the PD/LR and PD/PU modes must be done separately (and
  # for each of B0 and B1).  Note we only look for settled pointing (PO)
  # files.
  my @pcwpoFileList  = $filename->get('event', 'xrt', 'pcw?po', '*');
  my @wtwpoFileList  = $filename->get('event', 'xrt', 'wtw?po', '*');
  my @lrb0poFileList = $filename->get('event', 'xrt', 'lrb0po', '*');
  my @lrb1poFileList = $filename->get('event', 'xrt', 'lrb1po', '*');
  my @pub0poFileList = $filename->get('event', 'xrt', 'pub0po', '*');
  my @pub1poFileList = $filename->get('event', 'xrt', 'pub1po', '*');

  # get the output stem: sw or st and sequence number
  my $outstem = undef;
  foreach my $evt ( @wtwpoFileList,  @pcwpoFileList,
		    @lrb0poFileList, @lrb1poFileList,
		    @pub0poFileList, @pub1poFileList ) {
    $evt =~ /(s[wt]\d{11})x(pc|wt|lr|pu)/;
    if ( $1 ) {
      $outstem = $1;
      last;
    }
  }
  if (!defined $outstem) {
    $log->error(1, "failed to determine outstem to run XrtGrbLc, exit 1");
    return;
  }

  my $evtfiles = join(',', @pcwpoFileList, @wtwpoFileList); # all PC, WT modes
  $log->entry("XRT PC & WT modes event files found: $evtfiles");

  my $evtfiles_lrb0 = join(',', @lrb0poFileList); # PD/LR mode, B0
  $log->entry("XRT PD/LR mode B0 event files found: $evtfiles_lrb0");

  my $evtfiles_lrb1 = join(',', @lrb1poFileList); # PD/LR mode, B1
  $log->entry("XRT PD/LR mode B1 event files found: $evtfiles_lrb1");

  my $evtfiles_pub0 = join(',', @pub0poFileList); # PD/PU mode, B0
  $log->entry("XRT PD/PU mode B0 event files found: $evtfiles_pub0");

  my $evtfiles_pub1 = join(',', @pub1poFileList); # PD/PU mode, B1
  $log->entry("XRT PD/PU mode B1 event files found: $evtfiles_pub1");

  if ( (! $evtfiles) &&
       (! $evtfiles_lrb0) && (! $evtfiles_lrb1) &&
       (! $evtfiles_pub0) && (! $evtfiles_pub1) ) {
      $log->error(1, "No event files exist, exit 1");
      return ;
  }


  # xrthd, mkf files are used for PC, WT modes but not PD modes, so
  # shouldn't return if they're not present as we used to.
  my $xrthd_header = $filename->get('hk', 'xrt', 'hd');
  my $have_xrthd = -e $xrthd_header;
  if (! $have_xrthd ) {
    $log->error(1, "$xrthd_header xrthd file does not exist, can't process XRT PC, WT modes");
  }

  my $mkf_filter = $filename->get('filter', 's');
  my $have_mkf = -e $mkf_filter;
  if (! $have_mkf ) {
    $log->error(1, "$mkf_filter mkf file does not exist, can't process XRT PC, WT modes");
  }


  # Attitude file: pat if available, else sat.
  # Return if none found, but it's not used for the PD modes so maybe we
  # shouldn't?
  my $attitude = $filename->get('attcorr', 'p');
  if (-e $attitude) {
      $log->entry("Got pat attitude file $attitude");
  } else {
      $log->error(1, "$attitude pat attitude file does not exist, trying sat file");
      $attitude = $filename->get('attitude', 's');
      if (-e $attitude) {
	  $log->entry("Got sat attitude file $attitude");
      } else {
	  $log->error(1, "Unable to find attitude file, exit 1");
	  return ;
      }
  }

  $self->{XRT_HDFILE} = $xrthd_header;
  $self->{XRT_ATTFILE} = $attitude;

################################
# ADD KEYWORDS
################################
## (why do we have to go through all this to make the list of files??)
#   my @evtfiles1  = split(/","/, "$evtfiles");
#   my @evtfilesLR = split(/","/, "$evtfiles_lr");
#   my @evtfilesPU = split(/","/, "$evtfiles_pu");
#   my $evtfiles2  = join(',', @evtfiles1, @evtfilesLR, @evtfilesPU, $attitude);
#
#   my @evtfiles3 = split(/,/, $evtfiles2);
#
#   $self->write_standard_keyword(@evtfiles3);

  $self->write_standard_keyword( ( @wtwpoFileList,  @pcwpoFileList,
				   @lrb0poFileList, @lrb1poFileList,
				   @pub0poFileList, @pub1poFileList,
				   $attitude ) );

#################################ADD KEYWORDS


  ########################
  # GET GRB COORDINATES
  ########################

  # Read coords from job.par, but note these are never actually used:
  # see comments below on choosing "best" coords.
  my $ranom   =  $jobpar->read("burst_ra");
  my $decnom  =  $jobpar->read("burst_dec");

  # First try Lorella's catalog; if it fails,
  # then try JD's catalog; then try local catalog;
  # and finally send an email.

  # Read local catalog
  my $refPL = $self->Subs::UvotProduct::get_grb_coord_local();

  # Check if this is a GRB, if not ($refPL == -1), simply return
  # (Message is logged by get_grb_coord_local.)
  if (defined $refPL and $refPL == -1){
    return;
  }

  # Lorella's catalog.
  my $refP = $self->Subs::UvotProduct::get_swiftgrb_coord();
  my $catRead = "Lorella";

  # JD's catalog
  if (!defined $refP) {
      $refP = $self->Subs::UvotProduct::get_grb_coord();
      $catRead = "JD";
  }


  if (!defined $refP and !defined $refPL) {

    my $enddate = $jobpar->read("enddate");
    my $endtime = $jobpar->read("endtime");

    my $dateobj = Util::Date->new($enddate, $endtime);
    my $emjd = $dateobj->mjd();

    my $ndateobj = Util::Date->new();
    my $nmjd = $ndateobj->mjd();

    my $alerT    = $procpar->read("alerTime");

    if (!defined $alerT ) {
      $alerT = 2.0;
    }

    my $dt = $nmjd - $emjd;


    if ($dt > $alerT) {
	my $object = $jobpar->read('object');
	my $trigid = $jobpar->read('sequence');
	my $target = $jobpar->read('target');
	my $seqprocnum = $jobpar->read('seqprocnum');
	my $job_title  = $jobpar->read('job_title');
	my $obsdate    = $jobpar->read('obsdate');

	my $subject = "GRB not found in catalogs: $trigid";

	my $msg = "Unable to find GRB info in all catalogs for:\n\n" .
	    "object= $object\n" .
	    "sequence= $trigid\n" .
	    "target= $target\n" .
	    "processing number= $seqprocnum\n" .
	    "observation date= $obsdate\n" .
	    "job title= $job_title\n\n" .
	    "\tPlease, make sure that such info is actually available.\n\n" .
	    "\tYou can add the required info to the local GRB catalog\n" .
	    "by using\n\n" .
	    "/data/sdc/local/data/sdc4/apsop/Processman/PMT/bin/AddGRB.pl\n\n" .
	    "(Email generated in XrtGrbLc::make_xrt_grb_lc)\n";

	$self->sendEmail($subject, $msg, $watchers);
	$log->error(1, "Unable to find GRB in catalogs. Email notice sent to $watchers");

    } else {
	$log->error(1, "Unable to find GRB in catalogs. Email notice not sent, too soon");
    }

    return;
  }


  my @rItems = qw/bat_ra bat_dec bat_pos_err xrt_ra xrt_dec xrt_pos_err
      uvot_ra uvot_dec uvot_pos_err ot_ra ot_dec ot_pos_err/;

  if (defined $refPL) {

      if (!defined $refP) {
	  foreach my $it (@rItems) {
	      $refP->{$it} = $refPL->{$it};
	  }
	  $log->entry("Taking all coordinates from Local Catalog");
      } else {
	  if ((defined $refPL->{bat_over} and $refPL->{bat_over} == 1) or
	      (!defined $refP->{bat_ra} or !defined $refP->{bat_dec})) {
	      $refP->{bat_ra}      = $refPL->{bat_ra};
	      $refP->{bat_dec}     = $refPL->{bat_dec};
	      $refP->{bat_pos_err} = $refPL->{bat_pos_err};
	      $log->entry("BAT invalid in $catRead cat, or Local override");
	      $log->entry("Using BAT coordinates from Local Catalog");
	  } else {
	      $log->entry("Using BAT coordinates from $catRead Catalog");
	  }

	  if ((defined $refPL->{xrt_over} and $refPL->{xrt_over} == 1) or
	      (!defined $refP->{xrt_ra} or !defined $refP->{xrt_dec})) {
	      $refP->{xrt_ra}      = $refPL->{xrt_ra};
	      $refP->{xrt_dec}     = $refPL->{xrt_dec};
	      $refP->{xrt_pos_err} = $refPL->{xrt_pos_err};
	      $log->entry("XRT invalid in $catRead cat, or Local override");
	      $log->entry("Using XRT coordinates from Local Catalog");
	  } else {
	      $log->entry("Using XRT coordinates from $catRead Catalog");
	  }

	  if ((defined $refPL->{ot_over} and $refPL->{ot_over} == 1) or
	      (!defined $refP->{ot_ra} or !defined $refP->{ot_dec})) {
	      $refP->{ot_ra}      = $refPL->{ot_ra};
	      $refP->{ot_dec}     = $refPL->{ot_dec};
	      $refP->{ot_pos_err} = $refPL->{ot_pos_err};
	      $log->entry("OPT invalid in $catRead cat, or Local override");
	      $log->entry("Using OPT coordinates from Local Catalog");
	  } else {
	      $log->entry("Using OPT coordinates from $catRead Catalog");
	  }

	  if ((defined $refPL->{uvot_over} and $refPL->{uvot_over} == 1) or
	      (!defined $refP->{uvot_ra} or !defined $refP->{uvot_dec})) {
	      $refP->{uvot_ra}      = $refPL->{uvot_ra};
	      $refP->{uvot_dec}     = $refPL->{uvot_dec};
	      $refP->{uvot_pos_err} = $refPL->{uvot_pos_err};
	      $log->entry("UVOT invalid in $catRead cat, or Local override");
	      $log->entry("Using UVOT coordinates from Local Catalog");
	  } else {
	      $log->entry("Using UVOT coordinates from $catRead Catalog");
	  }
      }
  } else {
      $log->entry("Not in Local Catalog.  Taking all coordinates from $catRead Catalog");
  }


  # The GRB position must be chosen based on the "best" coordinates. If UVOT
  # coordinate are defined they will be used, else the Optical or OT
  # coordinates will be used if defined (with a warning email sent if they
  # disagree too much with the XRT ones), and finally if the XRT coordinates
  # are defined they will be used.  No light curve is made if only BAT
  # coordinates are available.

  if (defined $refP->{uvot_ra} && defined $refP->{uvot_dec} &&
      (!defined $refPL or !defined $refPL->{uvot_bad_data} or
       $refPL->{uvot_bad_data} == 0)) {

      $ranom  = $refP->{uvot_ra};
      $decnom = $refP->{uvot_dec};
      $log->entry("Will be using UVOT coordinates ($ranom, $decnom)");

  } elsif (defined $refP->{ot_ra} && defined $refP->{ot_dec} &&
	   (!defined $refPL or !defined $refPL->{ot_bad_data} or
	    $refPL->{ot_bad_data} == 0)) {
      $ranom  = $refP->{ot_ra};
      $decnom = $refP->{ot_dec};
      $log->entry("Will be using OPTICAL coordinates ($ranom, $decnom)");

      if (defined $refP->{xrt_ra} && defined $refP->{xrt_dec} &&
	  defined $refP->{xrt_pos_err}) {

          my $mincos = cos_d(2.0*$refP->{xrt_pos_err}/3600.0);
          my $costheta = sin_d($decnom) * sin_d($refP->{xrt_dec}) +
              cos_d($decnom) * cos_d($refP->{xrt_dec}) * cos_d($ranom - $refP->{xrt_ra});

	  if ($costheta < $mincos) {

	      my $subject = "OPTICAL GRB coordinates outside XRT coordinates";

	      my $msg = "Object = $object\nSequence = $seq\n\n"
		  . "Warning: While using OPTICAL coordinates for GRB, it appears that\n"
		  . "these coordinates ($ranom, $decnom) are not in agreement\n"
		  . "with the XRT ones ("
		  . $refP->{xrt_ra} . ',' . $refP->{xrt_dec} . ") within 2.*xrt_pos_err of "
		  . $refP->{xrt_pos_err} . " arcsec.\n"
		  . "\n\tPlease, verify that the catalog values are as expected.\n"
		  . "\tPositions from $catRead catalog unless Local overrides.\n";

	      $self->sendEmail($subject, $msg, $watchers);

	      $log->entry("make_xrt_grb_lc: Warning sent for Optical/XRT disagreement > 2*xrt_pos_err\n"
			  . "   xrt_ra=" . $refP->{xrt_ra} . " deg,  xrt_dec=" . $refP->{xrt_dec}
			  . " deg,  xrt_pos_err=" . $refP->{xrt_pos_err} . " arcsec" );
	  }
      }

  } elsif (defined $refP->{xrt_ra} && defined $refP->{xrt_dec} &&
	   (!defined $refPL or !defined $refPL->{xrt_bad_data} or
	    $refPL->{xrt_bad_data} == 0)) {
      $ranom  = $refP->{xrt_ra};
      $decnom = $refP->{xrt_dec};
      $log->entry("Will be using XRT coordinates ($ranom, $decnom)");

  } else {
      $log->error(1, "Only BAT coordinates, or none, are available.  No XRT light curve for GRB will be made.");
      return;
  }


  my $plotfile0 = $filename->get("lcplot", "xrt", "sr", 0);
  my $plotfile1 = $filename->get("lcplot", "xrt", "sr", 1);
  my $plotfile2 = $filename->get("lcplot", "xrt", "sr", 2);


  #############################
  #
  # MAKE XRT GRB LIGHT CURVE
  # PC AND WT MODES
  #
  #############################

  if ( $evtfiles && $have_xrthd && $have_mkf ) {

      $log->entry('Running xrtgrblc for PC & WT modes');

      my $odir = './tmp_dir_pcwt';
      if (-e $odir) {
	  system("rm -fr $odir");
      }

      # xrtgrblc needs one of these for every evtfile, even if they're all
      # the same.
      my $num1 = $#pcwpoFileList + 1;
      my $num2 = $#wtwpoFileList + 1;
      my $mkffiles = join(',', ($mkf_filter)x($num1 + $num2));
      my $xhdfiles = join(',', ($xrthd_header)x($num1 + $num2));
      my $attfiles = join(',', ($attitude)x($num1 + $num2));

      my $xrtgrblc=Util::HEAdas->new("xrtgrblc")->is_script( 1 );
      $xrtgrblc->params({
	  evtfiles     => $evtfiles,
	  mkffiles     => $mkffiles,
	  xhdfiles     => $xhdfiles,
	  attfiles     => $attfiles,
	  outdir       => $odir,
	  outstem      => $outstem,
	  srcra        => $ranom,
	  srcdec       => $decnom,
	  minenergy    => 0.3,
	  maxenergy    => 10.0,
	  splitenergy  => 0.0,
	  splitorbits  => 'no',
	  cutlowbins   => 'yes',
	  plotftype    => '/gif',
	  clobber      => 'yes',
	  history      => 'yes',
	  cleanspec    => 'yes',
	  # 2015-02-24 REW added lccorravg, pcreglist, wtreglist per 2015-01-14 email from Matteo
	  lccorravg    => 'no',
	  pcreglist    => 2,
	  wtreglist    => 2,
	  });

      $xrtgrblc->run();


      # Replot the light curve. A better title should be added, and
      # the X-axis label should not say "BAT Trigger" if no BAT trigger

      # Setup some plot labels
      my $trigtime = $jobpar->read('burst_time');
      my $trigtimeSrc = "burst_time Parameter";
      my $trigdate = undef;
      my $xlab;

      # Note: trigger time read here appears to only be used for the axis label
      # of the xsr_lc plot below!
      if (!defined $trigtime or $trigtime == 0){
	  $trigtime = $self->Subs::UvotProduct::getTrigFromDB();
	  $trigtimeSrc = "TDRSS Catalog";

	  if (!defined $trigtime){
	      $trigtime = $self->Subs::UvotProduct::getTrigFromSWIFTCatalog();
	      $trigtimeSrc = "Lorella Catalog";

	      if (!defined $trigtime){
		  $trigtime = $self->Subs::UvotProduct::getTrigFromJDCatalog();
		  $trigtimeSrc = "JD Catalog";

		  if (!defined $trigtime){
		      $trigtime = $self->Subs::UvotProduct::getTrigFromLocalCatalog();
		      $trigtimeSrc = "Local Catalog";
		  }
	      }
	  }
      }


      if (defined $trigtime){
	  my $dateobj = Util::Date->new($trigtime);
	  $trigdate = $dateobj->date().'T'.$dateobj->time();
      }

      if (defined $trigtime and defined $trigdate){
	  $xlab = "Time(s) since BAT trigger time(UT $trigdate/MET $trigtime)";
      } else {
	  $xlab = "Time [MET(s)]";
      }


      my $evtfile = ( @pcwpoFileList, @wtwpoFileList )[ 0 ];
      my $evtfits = Util::FITSfile->new($evtfile);
      my $dateobs = $evtfits->keyword( 'DATE-OBS' );
      $dateobs =~ s/[T ]\d{2}:\d{2}:[\d\.]+$//;
      my $title = uc( $object ) . ' SWIFT XRT '
	  . $self->jobpar->read('sequence')
	  . ' (' . $dateobs . ')';

      # change to the dir where the qdp file is
      my $curdir  = getcwd( );
      chdir $odir;


      my $qdp = Util::Xanadu->new( 'qdp' );
      $qdp->verbose( 1 );
      $qdp->seriousness( 1 );

      $qdp->script(
		   "${outstem}_xpwetsrab.qdp",
		   '/null',
		   'SCR  WHITE',
		   'TIME OFF',
		   'LAB  F',
		   "LAB  T \"$title\"",
		   "LAB  X \"$xlab\"",
		   'YPL  OFF 1 3 4 5 7 8',
		   'LOG',
		   'GAP  ER',
		   'RES',
		   "cpd ${outstem}xsr_lc.gif/gif",
		   "PLOT",
		   "QUIT"
		   )->run( );


      chdir $curdir;

      #
      # plot the spectra
      #
      my @phaF = ();
      push @phaF, "$odir/${outstem}_xpcetsr.pha";
      push @phaF, "$odir/${outstem}_xwtetsr.pha";

      $self->plotSpectra( \@phaF, "${outstem}xsr_ph.gif/gif" ); # NOTE basename here

      #
      # move files into their final names
      # sr = source
      #
      rename "$odir/${outstem}_xpcetsr.pha", "$odir/${outstem}xpcsr.pha";
      rename "$odir/${outstem}_xpcet.arf",   "$odir/${outstem}xpcsr.arf";
      rename "$odir/${outstem}_xwtetsr.pha", "$odir/${outstem}xwtsr.pha";
      rename "$odir/${outstem}_xwtet.arf",   "$odir/${outstem}xwtsr.arf";
      rename "$odir/${outstem}_xpcetsra.lc", "$odir/${outstem}xpcsr.lc";
      rename "$odir/${outstem}_xwtetsra.lc", "$odir/${outstem}xwtsr.lc";

      # remove un-needed files
      unlink "$odir/${outstem}_xpcetbg.pha", "$odir/${outstem}_xwtetbg.pha",
      "$odir/${outstem}_xpwetsrfb_lc.gif", "$odir/${outstem}_xpwetsrcb_lc.gif",
      "$odir/${outstem}_info.fits";

      #unlink "$odir/${outstem}_xpwetsrab.qdp", "$odir/${outstem}_xpwetsrab.pco";

      # Fixup keywords: rewrite (or add) ANCRFILE with the arf file's new name,
      # on both HDUs.
      if ( -e "$odir/${outstem}xpcsr.pha" ) {
	  Util::FITSfile->new("$odir/${outstem}xpcsr.pha", 0)
	      ->keyword("ANCRFILE", "${outstem}xpcsr.arf");
	    Util::FITSfile->new("$odir/${outstem}xpcsr.pha", 1)
		->keyword("ANCRFILE", "${outstem}xpcsr.arf");
	}

      if ( -e "$odir/${outstem}xwtsr.pha" ) {
	  Util::FITSfile->new("$odir/${outstem}xwtsr.pha", 0)
	      ->keyword("ANCRFILE", "${outstem}xwtsr.arf");
	    Util::FITSfile->new("$odir/${outstem}xwtsr.pha", 1)
		->keyword("ANCRFILE", "${outstem}xwtsr.arf");
	}

      # Concat log files with the ones in the main directory
      $self->concatLogs( $outstem, $odir, $curdir );

      # Move any remaining files in $dir to $curdir.
      if (-e $odir) {
	  system("mv $odir/* $curdir");
	  system("rm -fr $odir");
      }

  } # end of PC, WT modes



  ##########################################################################
  #
  # MAKE XRT GRB LIGHT CURVES
  # PD MODES (LR and PU)
  #
  # Names for the output files try to follow the PC,WT convention above
  # (eg, "${outstem}xpcsr.lc") but need to include b0/b1.  PC,WT names
  # include "sr" ("Swift GRB Product File Naming Convention" says
  # sr=source); even though xrtproducts help says PD modes have no spatial
  # information and all events are used, Matteo Perri (XRT team at ASDC)
  # says "sr" is appropriate for these modes too because they were only
  # used for extremely bright sources.  "po" is not included because we're
  # using only the settled pointing event files.
  ##########################################################################

  if ( $evtfiles_lrb0 ) {
      $self->make_pd_mode_lc( $evtfiles_lrb0, 'xlrb0sr', $outstem );
  }

  if ( $evtfiles_lrb1 ) {
      $self->make_pd_mode_lc( $evtfiles_lrb1, 'xlrb1sr', $outstem );
  }

  if ( $evtfiles_pub0 ) {
      $self->make_pd_mode_lc( $evtfiles_pub0, 'xpub0sr', $outstem );
  }

  if ( $evtfiles_pub1 ) {
      $self->make_pd_mode_lc( $evtfiles_pub1, 'xpub1sr', $outstem );
  }


} # end of make_xrt_grb_lc method


##########################################################################
#
# make_pd_mode_lc( $evtfiles, $basename, $outstem )
#
# Run xrtproducts to make light curve, arf, and pha products for one of the
# XRT Photo Diode (PD/LR and PD/PU) modes: xlrb0, xlrb1, xpub0,
# xpub1. (b1=bias subtracted onboard, b0=not subtracted.)  They're all done
# the same way, but unlike PC and WT must be done separately; call this
# with the appropriate event file for each one.  Files created are *.lc,
# *.arf, *.pha, *_lc.gif, and *_ph.gif.  Work is done in a temporary
# directory, as for the PC, WT modes.  The procedure used is from Matteo
# Perri (ASDC) (emails 2014-03-28 and 2014-04-07).
#
# Parameters:
#   $evtfiles: String of input event file names for this mode, comma
#              separated if more than one.
#   $basename: Base name of the output files, will be appended onto $outstem.
#              Eg., 'xlrb0sr'.
#   $outstem:  Stem for the output file names.  Eg, 'sw00106415000'.
#
#   Eg., $basename='xlrb0sr', $outstem='sw00106415000' produces files named
#   sw00106415000xlrb0sr.lc, etc.
#
##########################################################################

sub make_pd_mode_lc {

    my $self = shift;
    my ($evtfiles, $basename, $outstem) = @_;
    my $log  = $self->log();

    my $attfile = $self->{XRT_ATTFILE};
    my $hdfile = $self->{XRT_HDFILE};
    if (not $hdfile or not -f "$hdfile") {
       $log->entry("Skipping xrtproducts for PD mode ${basename} without hdfile");
       return;
    }

    $log->entry("Running xrtproducts for PD mode ${basename}");

    my $curdir = getcwd( );
    my $odir = "./tmp_dir_${basename}";
    if (-e $odir) {
	system("rm -fr $odir");
    }

    # Names for the output files.

    my $lcName  = "${outstem}${basename}.lc";
    my $phaName = "${outstem}${basename}.pha";
    my $arfName = "${outstem}${basename}.arf";

    # Run xrtproducts to create the light curve and related files.  The
    # parameters listed were specified by Matteo Perri (mostly) as the
    # appropriate values for processing all the PD modes and should be the
    # complete set (email from Matteo 2014-04-07).

    # NOTE: Setting stemout=>"${outstem}${basename}" is a partial
    # workaround suggested by Matteo for a bug in xrtproducts v.0.4.1 and
    # earlier, which hardcodes the plot title in the light curve gif to
    # "Light curve (${stemout}pcsr.lc)" for ALL modes.  (The pha gif title
    # appears to use $phafile.)  This also changes the names of the .gif
    # files, but since we specify all the others explicitly, it appears to
    # have no other effect.  Change this back to stemout=>$outstem, and the
    # rename commands below, if Matteo gives us a patch for the bug.

    my $xrtproducts = Util::HEAdas->new("xrtproducts")->is_script( 1 );
    $xrtproducts->params( {
	infile     => $evtfiles,
	outdir     => $odir,
	lcfile     => $lcName,
	phafile    => $phaName,
	arffile    => $arfName,
	# This stemout partially works around the bug in lc.gif plot titles
	stemout    => "${outstem}${basename}",
	binsize    => -99,
	display    => 'no',
	gtifile    => 'NONE',
	rmffile    => 'CALDB',
	mirfile    => 'CALDB',
	transmfile => 'CALDB',
	inarfile   => 'CALDB',
	psfile     => 'CALDB',
	vigfile    => 'CALDB',
	psfflag    => 'yes',
	extended   => 'no',
	expofile   => 'NONE', # only applicable for PC/WT modes
        correctlc  => 'no',
	hdfile     => $hdfile,
	attfile    => $attfile,
	plotdevice => 'gif',
	pilow      =>   30,
	pihigh     => 1000,
	clobber    => 'yes',
	chatter    =>    3,
	history    => 'yes',
	cleanup    => 'yes'
	} );
    $xrtproducts->run();

    # Move files into their final names
    # Doesn't look like xrtproducts has a way to set the .gif names;
    # xrtproducts names the gif files ${stemout}lc.gif and ${stemout}ph.gif
    rename "${odir}/${outstem}${basename}lc.gif", "${odir}/${outstem}${basename}_lc.gif";
    rename "${odir}/${outstem}${basename}ph.gif", "${odir}/${outstem}${basename}_ph.gif";

    # remove unneeded files

    # Fixup ANCRFILE keywords in the .pha file with the .arf file's name
    if ( -e "${odir}/${phaName}" ) {
	my $phafile = Util::FITSfile->new("${odir}/${phaName}");
	$phafile->ext(0);
	$phafile->keyword("ANCRFILE", $arfName);
	$phafile->ext(1);
	$phafile->keyword("ANCRFILE", $arfName);
    }

    # Concat log files with the ones in the main directory
    $self->concatLogs( $outstem, $odir, $curdir );

    # Move any remaining files in $dir to $curdir.

    if (-e $odir) {
	system("mv $odir/* $curdir");
	system("rm -fr $odir");
    }

}  # end of make_pd_mode_lc method


###############################################################
#
# write_standard_keyword: Add standard keywords to FITS files.
# This one is currently called by XrtGrbLc and UvotProducts.
# There is a very similar routine in SW0WrapUp (this one was probably
# copied from there).
#
# Parameter:
#    @file1: List of filenames to add keywords to.
#
###############################################################

sub write_standard_keyword {

    my $self=shift;
    my $log      = $self->log();
    my $filename = $self->filename();
    my $jobpar   = $self->jobpar();
    my $procpar  = $self->procpar();
    my @file1 = @_;

    #######################
    # Setup for UTCF value
    #######################
    my $time_file = $filename->get('timedata', 'swift', '', 0);
    my (@utcf_times, @utcf);
    if( -f $time_file ){
      my $time_fits = Util::FITSfile->new($time_file, 'UTCF');
      @utcf_times = $time_fits->cols('TIME')->table();
      @utcf = $time_fits->cols('UTCF')->table();
    }

    my $soft_version = $jobpar->read('softver');
    my $cal_version = $jobpar->read('caldbver');
    my $reprocess = $jobpar->read('reprocess');

    ############################
    # Trigger time, if relevant
    ############################
    my $trigtime = $jobpar->read('burst_time');


    ####################################
    # loop over all FITS files
    ####################################
    my $file;

    foreach $file (@file1 ) {
#        $log->entry("Doing file $file.");
	$log->entry("Adding standard keywords to file $file");

        my $fits = Util::FITSfile->new($file);
	my $is_tdrss = $file =~ /s[wt][t\d]\d{10}ms/;
	my $is_batevt = $file =~ /s[wt][t\d]\d{10}bev/;
	my $is_eng = $file =~ /s[wt][t\d]\d{10}.*en\.hk$/;
	my $is_att = $file =~ /s[wt][t\d]\d{10}.at\.fits$/;
	my $is_sm = $file =~ /s[wt][t\d]\d{10}bsmcb\.fits$/;

        ##################################
        # loop over HDUs in the FITS file
        ##################################
        my $nhdus = $fits->nhdus();
	unless ($nhdus) {
#	  $log->error(2, "Cannot get number of HDUs for FITS file $file.");
	  next;
	}

        my $hdu;
        for($hdu=0; $hdu<$nhdus; $hdu++) {
            $fits->ext($hdu);
	    my $extname = $hdu==0 ? '' : $fits->keyword('EXTNAME');

            ################################
            # write keywords to the file
            ################################
            $fits->begin_many_keywords();

            $fits->keyword('PROCVER', $procpar->read('version'),
                           'Processing script version' );

	    unless( $is_eng || ($extname && $extname=~/GTI/) ){
	      $fits->keyword('SOFTVER', $soft_version, 'HEASOFT and Swift versions');
	      $fits->keyword('CALDBVER', $cal_version, 'CALDB index versions used');
	    }

	    if($hdu==0){
	      $fits->keyword('TIMESYS', 'TT', 'time system');
	      $fits->keyword('MJDREFI', 51910, 'MJD reference day 01 Jan 2001 00:00:00');
	      $fits->keyword('MJDREFF', 7.428703700000000E-04,
			     'MJD reference (fraction of day) 01 Jan 2001 00:00:00');
	      $fits->keyword('CLOCKAPP', 'F', 'If clock correction are applied (F/T)');
	      $fits->keyword('TIMEUNIT', 's', 'Time unit for timing header keywords');

	      $fits->keyword('REPROC', 'T', 'Is this from a bulk reprocessing run?') if $reprocess eq 'yes';

	    }else{
	      $fits->keyword('TIERRELA', '1.0E-8', '[s/s] relative errors expressed as rate');
	      $fits->keyword('TIERABSO', '1.0', '[s] timing  precision  in seconds');
	    }

	    unless($is_tdrss || $is_sm){
              $fits->keyword('OBS_ID', sprintf("'%011d'", $jobpar->read('sequence')),
                             'Observation ID' );

              $fits->keyword('SEQPNUM', int($jobpar->read('seqprocnum')),
                             'Number of times the dataset processed' );
	      $fits->keyword('TARG_ID', int(sprintf("%d",$jobpar->read('target'))),
			     'Target ID');
	      $fits->keyword('SEG_NUM', int(sprintf("%d",$jobpar->read('obs'))),
			     'Segment number');

	      $fits->keyword('OBJECT', $jobpar->read('object'), 'Object name');
	      $fits->keyword('RA_OBJ', $jobpar->read('burst_ra'), '[deg] R.A. Object');
	      $fits->keyword('DEC_OBJ', $jobpar->read('burst_dec'), '[deg] Dec Object');

	      $fits->keyword('RA_PNT', $jobpar->read('ra'), '[deg] RA pointing');
	      $fits->keyword('DEC_PNT', $jobpar->read('dec'), '[deg] Dec pointing');
	      $fits->keyword('PA_PNT', $jobpar->read('roll'), '[deg] Position angle (roll)');

	      $fits->keyword('TRIGTIME', $trigtime, '[s] MET TRIGger Time for Automatic Target')
		if $trigtime;
	    }

	    if($is_batevt){
	      $fits->keyword('CATSRC',
			     $jobpar->read('burst_cat_src') eq 'yes' ? 'T' : 'F');
	    }

	    if ( $is_att ) {
	      $fits->keyword('-ATTSTATU', ' ');
	    } else {

	      #############################
	      # Determine and set UTCFINIT
	      #############################
	      my $tstart = $fits->keyword('TSTART');
	      unless( $tstart ){
		my $date = $fits->keyword('DATE-OBS');
		if( $date ){
		  my $start = Util::Date->new($date);
		  $tstart = $start->seconds();
		}
	      }

	      if( $tstart && @utcf ){
		my $itime = 0;
		while( $tstart > $utcf_times[$itime] &&
		       $itime < $#utcf_times){ $itime++ }
		$fits->keyword('UTCFINIT', $utcf[$itime], '[s] UTCF at TSTART');
	      }
	    }

            $fits->end_many_keywords();
	  }

      }

} # end of write_standard_keyword method


###############################################################
#
# sendEmail: Send an email
#
# Parameters:
#    $subject:  Subject
#    $msg:      Text of the message, can contain \n etc.
#    $watchers: Recipient(s), email addresses separated by spaces
#
# Apostrophes will be removed from subject, because they break the
# shell command line we construct.  They're OK within the msg text.
###############################################################

sub sendEmail {

  my $self=shift;
  my ($subject, $msg, $watchers) = @_;


  my $jobpar   = $self->jobpar();

  my $trigid =$jobpar->read('sequence');

  my $tm = time();
  my $tfile = "/tmp/msg_$tm".'_'.$trigid;

  # Remove any apostrophes from subject, because
  # they will break the constructed command line.
  $subject =~ s/\'//g ;

  if (open OUTF, ">$tfile") {
    print OUTF "$msg";
    close OUTF;

    my $hostname = `hostname`;
    chomp $hostname;

    if ($hostname =~ /^sdc/ ) {

      my $cmd = "mail -s \'$subject\' $watchers < $tfile";
      `$cmd`;
      unlink $tfile;

    } else {

      my $rv = system("scp -q $tfile apsop\@sdc:/tmp  >& /dev/null");

      if ($rv == 0) {

	my $cmd1 = "ssh -nq apsop\@sdc \"mail -s \'$subject\' $watchers < $tfile\"";
	my $rval = system($cmd1);

	unlink $tfile;
	system("ssh -nq apsop\@sdc rm -f $tfile");
      }
    }
  }
}  # end of sendEmail method



###########################################################
#
# plotSpectra - groups and plots spectra.
# Runs grppha and xspec.
#
# Parameters:
#    $refspectra: Reference to array of .pha files to plot.
#    $outplot:    Name of the resulting plot file.  Deleted if
#                 no plot is made (ie, plot file comes out empty).
#
############################################################

sub plotSpectra {

  my $self = shift;

  my ( $refspectra, $outplot ) = @_;

  my $log      = $self->log();

  my $jobpar   = $self->jobpar();

  my $procpar  = $self->procpar();

  my $watchers = $procpar->read("watchers");
  my $headas   = $procpar->read('headas');

  my $GRB      =  $jobpar->read('object');
  $GRB         =  uc($GRB);

  my @localclean = ( );   # list of files we'll delete before returning

  # check which spectra we actually have
  my @spectra = ( );
  foreach my $spec1 ( @$refspectra ) {
    if ( -e $spec1 ) {
      push @spectra, basename( $spec1 );
    }
  }

  return unless @spectra;


  # change to the dir where the spectra are
  my $workdir = dirname( $refspectra->[ 0 ] );
  my $curdir  = getcwd( );
  chdir $workdir;

  #
  # Get the response matrices and group each spectrum
  #
  my $i = 1;
  my @grpspectra = ( );
  my $grppha=Util::HEAdas->new("grppha");
  my @color = ();

  foreach my $spec ( @spectra ) {
      # Get name of response matrix file (RMF) from RESPFILE keyword
      # (FITSfile->keyword strips quotes and whitespace (inc. \n))
      # If can't, throw error and go on to next spectrum.
      my $resp = Util::FITSfile->new($spec,1)->keyword('RESPFILE');
      if ( $resp ) {    # defined and not null
	  $log->entry("XrtGrbLc plotSpectra: $spec RESPFILE is $resp\n");
      } else {
	  $log->error(2, "Could not read keyword RESPFILE from $spec");
	  next;
      }

      # Copy RMF file from caldb.
      # Error and go on to next if we can't (there may be more info in the
      # proclog file).
      # Should be more general than this, but RMF should be in current
      # directory.
      my $caldbResp = "$ENV{CALDB}/data/swift/xrt/cpf/rmf/$resp";
      if ( -e $caldbResp ) {
	  my $trv = system( "cp $caldbResp ./$resp" );
	  if ($trv == 0) {
	      $log->entry("got RMF file $resp from CALDB.\n");
	      push @localclean, $resp;
	  } else {
	      $log->error(2,"Could not copy $caldbResp to $workdir, status=$trv");
	      next;
	  }
      } else {
	  $log->error(2,"RMF file $caldbResp not found in CALDB!");
	  next;
      }

      # Run grppha to group
      my $grpspec = $spec;
      $grpspec =~ s/\.(pha|pi)$/_grp.$1/;

      $grppha->params( {
	  infile   => $spec,
	  outfile  => $grpspec,
	  comm     => 'group min 20',
	  tempc    => 'exit',
	  clobber  => 'yes'
	  } );

      $grppha->run();


      push @grpspectra, $grpspec;
      push @localclean, $grpspec;

      # Build up color command for plotting
      $grpspec =~ /x(pc|wt)/;
      my $mode = $1;
      push @color, sprintf( "setplot command color %d on %d",
			    $mode eq 'pc' ? 2 : 4, $i );

      $i++;
  }

  #
  # Write and run the xspec script
  # A prettier plot should be created (using setplot command)
  #

  # Setup some plot labels
  my $fits = Util::FITSfile->new($grpspectra[ 0 ]);
  my $dateobs = $fits->keyword( 'DATE-OBS' );
  $dateobs =~ s/[T ]\d{2}:\d{2}:[\d\.]+$//;
  my $title = $GRB . ' SWIFT XRT '
            . $self->jobpar->read('sequence')
            . ' (' . $dateobs . ')';

  my $xspec = Util::Xanadu->new("xspec");

  $xspec->verbose( 1 );
  $xspec->seriousness( 1 );

  # Build the xspec command list
  my @commands = ();

  push @commands, ('query no',
                   'setplot splashpage off',
                   '/*',
		   'setplot energy',
                   'setplot command LA OT " "',
		   "setplot command LA T $title",
		   'setplot command TIME OFF',
		   'setplot command SCR  WHITE',
		   'setplot command GAP  ER',
		   'setplot command RES',
		  );

  foreach my $c (@color) {
      if ($c !~ /^\s*$/) {
	  push @commands, $c;
      }
  }

  $i = 1;
  foreach my $spec ( @grpspectra ) {
      push @commands, ("data $i:$i $spec",
		       '/*',
		       '/*',
		       'ignore bad',
		       "ignore $i:**-0.3 10.0-**",
		       );
      $i++;
  }

  push @commands, ("cpd $outplot",
		   'plot ldata',
		   'exit'
		   );

  $log->entry("XrtGrbLc plotSpectra:  xspec command list:\n");
  $log->text( join("\n", @commands) );

  # Execute xspec
  $xspec->script(@commands)->run();


  my $gif = (split/\//, $outplot)[0];
  if (!-s $gif) {
      $log->error(1, "Attempted to create plot file $gif but its size is 0. File removed");
      unlink $gif;
  }

  # Remove files no longer needed and return to main directory
  unlink @localclean;
  chdir $curdir;

}  # end of plotSpectra method

##################################################################
#
# concatLogs( $outstem, $odir, $curdir )
#
# Concatenate the log files in $odir (./tmp_dir_*) to the ones in the main
# working directory $curdir.  (Output buffering should not be a problem,
# because Log.pm opens and closes the files for each entry written.)  The
# file names all begin with $outstem.  Removes the log files from $odir.
#
################################################################
sub concatLogs {

    my $self = shift;
    my ($outstem, $odir, $curdir) = @_;

    if (-e $odir) {

	if (-e "$odir/${outstem}pjl.html"){
	    `cat $odir/${outstem}pjl.html >> $curdir/${outstem}pjl.html`;
	    unlink "$odir/${outstem}pjl.html";
	}

	if (-e "$odir/${outstem}per.html"){
	    `cat $odir/${outstem}per.html >> $curdir/${outstem}per.html`;
	    unlink "$odir/${outstem}per.html";
	}

	if (-e "$odir/${outstem}pin.html"){
	    `cat $odir/${outstem}pin.html >> $curdir/${outstem}pin.html`;
	    unlink "$odir/${outstem}pin.html";
	}

    }

    return;
}  # end of concatLogs method


sub deg2rad
{
    my ($deg) = @_;
    return PI * $deg / 180;
}

sub sin_d
{
    my ($deg) = @_;
    return sin(deg2rad($deg));
}

sub cos_d
{
    my ($deg) = @_;
    return cos(deg2rad($deg));
}


1;