Subs::BATLightCurve (version $)


package Subs::BATLightCurve;
##############################################################################
#
# DESCRIPTION: Produce BAT lightcurve gif file
#
# HISTORY: $Log: BATLightCurve.pm,v $
# HISTORY: Revision 1.8  2014/08/14 09:23:45  apsop
# HISTORY: VERSION comment and minor formatting only (other changes were reverted).
# HISTORY:
# HISTORY: Revision 1.7  2012/08/15 07:56:28  apsop
# HISTORY: Discard output gif if it's empty. [Bob W.]
# HISTORY:
# HISTORY: Revision 1.6  2012/01/12 06:52:03  apsop
# HISTORY: Changes going to proc3.15.03
# HISTORY:
#
# VERSION: $Revision: 1.8 $
#
##############################################################################

use Subs::Sub;

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

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

        $self->{DESCRIPTION}= 'Make BAT light curves gif file';
        return $self;
}

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

sub body {
	my $self = shift;

    my $log      = $self->log;
    my $filename = $self->filename;
    $self->gif_light_curve;
}

sub gif_light_curve {
  my $self = shift;


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


  my $title = uc($jobpar->read('object')) .' SWIFT BAT '.
    $jobpar->read('sequence') .' ('.$jobpar->read('obsdate').')';

  my $log      = $self->log;
  my $filename = $self->filename;
  $log->entry('Making BAT light curves gif file');

	######################################
	#  check if files we needed exist
	######################################

  my $infile1 = $filename->get('lightcurve', 'bat', 'ev1s', 0);
  if (! -e $infile1 ) {
    $log->error(1, "$infile1 file does not exist, bailing");
    return ;
  }

  my $infile2 = $filename->get('lightcurve', 'bat', 'evms', 0);
  if (! -e $infile2 ) {
    $log->error(1, "$infile2 file does not exist, bailing");
    return ;
  }

  my $msbce = $filename->get('tdmess', 'bat', 'ce', 0);
  if (! -e $msbce ) {
    $log->error(1, "$msbce file does not exist, no DATETRIG,TRIGTIME,OBJECT");
  }

  my $bevbu = $filename->get('burstgti', 'bat', '', 0);
  if (! -e $bevbu ) {
    $log->error(1, "$bevbu file does not exist, no T90DUR,TOTSTART,TOTSTOP");
  }

  my $tmp = substr($infile1, 0, 18);
  my $out_qdp = $tmp."bevms.qdp";
  my $out_gif=$filename->corresponding("lightcurve", "lcplot", $infile1);

  $log->entry("BAT lightcurve input=$infile1, $infile2 output= $out_gif");

  my $bevbu1 ="";
  my $msbce1   = Util::FITSfile->new($msbce,0);
  my ($trigdate,$trigtime,$t90dur,$totstart,$totstop,$ext,$object)=0;
  $trigdate = $msbce1->keyword('DATETRIG');
  $trigtime = $msbce1->keyword('TRIGTIME');
  $object   = $msbce1->keyword('OBJECT');
  for ( $ext=0 ;$ext<11 ; $ext++ ) {
    $bevbu1 = Util::FITSfile->new($bevbu,$ext);
    $t90dur   = $bevbu1->keyword('T90DUR');
    $totstart = $bevbu1->keyword('TOTSTART');
    $totstop  = $bevbu1->keyword('TOTSTOP');
    if ($t90dur && $totstart && $totstop ) {
      last;
    }
  }
  if (!$object) {
    $object="GRB";
  }

  my ($t1,$t2,$t3)=0;
  if (defined($totstart) && defined($totstop) &&
      defined($t90dur)   && defined($trigtime)) {
    $t1=$totstart-0.5*$t90dur-$trigtime;
    $t2=$totstop +0.5*$t90dur-$trigtime;
    if ($t90dur < 1.5 ) {
      $t1=$totstart-$trigtime-1.5;
      $t2=$totstart-$trigtime+1.5;
    }
    $t3=$t1+0.8*($t2-$t1);
  } else {
    $t1=-5;
    $t2=+5;
    $t3=$t1+0.8*($t2-$t1);
  }

  ###########################
  #  make qdp file
  ###########################

  my $fits1 = Util::FITSfile->new($infile1, 1); #ext 1 has lightcurve
  my @time1 = $fits1->cols("TIME")->table();
  my @rate1 = $fits1->cols("RATE")->table();
  my @err1 = $fits1->cols("ERROR")->table();

  my $fits2 = Util::FITSfile->new($infile2, 1); #ext 1 has lightcurve
  my @time2 = $fits2->cols("TIME")->table();
  my @rate2 = $fits2->cols("RATE")->table();
  my @err2 = $fits2->cols("ERROR")->table();

  open FILE, ">$out_qdp" or die "failed to open $out_qdp for writing\n";

  print  FILE "read serr 2\n";
  print  FILE "skip single\n";

  # initialize variables in list to zero
  $_ = 0 for my ($t,$err,$r,$max1,$max2,$min1,$min2,$maxerr,$minerr);
  my $i=0; $max1=$min1=$rate1[0];
  foreach $t (@time1) {
    $err=sqrt($err1[$i]**2+$err1[$i+1]**2+$err1[$i+2]**2+$err1[$i+3]**2);
    $r=$rate1[$i]+$rate1[$i+1]+$rate1[$i+2]+$rate1[$i+3];
    if ($max1 < $r) {
      $max1=$r; $maxerr=$err;
    }
    if ($min1 > $r) {
      $min1=$r; $minerr=$err;
    }
    printf FILE "%13.6G %13.6G %13.6G\n",$t-$trigtime, $r, $err;
    $i=$i+4;
  }
  $min1=$min1-$minerr;
  $max1=$max1+$maxerr;

  print FILE "NO NO NO\n";

  $i=0; $max2=$min2=$rate2[0];
  foreach $t (@time2) {
    $err=sqrt($err2[$i]**2+$err2[$i+1]**2+$err2[$i+2]**2+$err2[$i+3]**2);
    $r=$rate2[$i]+$rate2[$i+1]+$rate2[$i+2]+$rate2[$i+3];
    if ($max2 < $r) {
      $max2=$r; $maxerr=$err;
    }
    if ($min2 > $r) {
      $min2=$r; $minerr=$err;
    }
    printf FILE "%13.6G %13.6G %13.6G\n",$t-$trigtime, $r, $err;
    $i=$i+4;
  }
  $min2=$min2-$minerr;
  $max2=$max2+$maxerr;

  close FILE;
  my $max11=$max1*0.9;
  my $max22=$max2*0.9;
  my $max111=$max1*1.1;
  my $max222=$max2*1.1;
  my $min111=$min1*1.1;
  my $min222=$min2*1.1;

  ##################################
  #  plot lightcurve to gif file
  ##################################

  ################################################
  # Setup a Xanadu QDP object
  ################################################
  my $qdp = Util::Xanadu->new( 'qdp' );
  $qdp->verbose( 0 );
  $qdp->seriousness( 1 );
  ###################################################
  # Write the QDP script - the first line MUST be
  # the name of the '.qdp' file to read from and the
  # second line MUST be the initial plot device
  # (should be '/null')
  ####################################################
  $qdp->script(
	       "$out_qdp",
	       "/null",
	       "RESCALE X $t1,$t2",
	       "WIN    1",
	       "RESCALE Y $min111 $max111",
	       "LINE ON ON 1",
	       "LABEL 1 POS $t3 $max11 \"15-350 keV - 1 second\"",
	       "LABEL 2 POS -$t1 0",
	       "LABEL 2 LI 0 10 LS 2",
#	       "LAB  T  Swift-BAT $object",
	       "LAB  T  $title",
	       "LAB  Y  Counts/sec/det",
	       "WIN      2" ,
	       "RESCALE Y $min222 $max222",
	       "LINE OFF ON 2",
	       "LABEL 3 CO 2 POS $t3 $max22 \"15-350 keV - 64 msec\"",
	       "LABEL 4 POS -$t1 0",
	       "LABEL 4 LI 0 10 LS 2",
	       "LAB X Time(s) since BAT trigger time(UT $trigdate/MET $trigtime)",
	       "LAB  Y  Counts/sec/det",
	       "WIN  ALL",
	       "SCR       WHITE",
	       "LABEL     FILE",
	       "TIME      OFF",
	       "PLOT      VERTICAL",
	       "CPD      $out_gif/gif",
	       "PLOT",
	       "QUIT"
	      )->run( );

  if ( -e $out_qdp) {
    unlink $out_qdp;
  }

  if (-z $out_gif ) {
    $log->error(1, "$out_gif is empty, discarding");
    unlink($out_gif);
  }

}				#end of sub gif_light_curve