Subs::BATLightCurve (version 0.0)

package Subs::BATLightCurve;
##############################################################################
#
# DESCRIPTION: Produce BAT lightcurve gif file
#
#
# VERSION: 0.0
#
##############################################################################

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 $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 not exist, exit 1");
	 	return ;
     }
     my $infile2 = $filename->get('lightcurve', 'bat', 'evms', 0);
     if(! -e $infile2 ){
       	$log->error(1, "$infile2 file not exist, exit 1");
	 	return ;
     }
	my $msbce = $filename->get('tdmess', 'bat', 'ce', 0);
    if(! -e $msbce ){
    	$log->error(1, "$msbce file not exist, no DATETRIG,TRIGTIME,OBJECT");
	}
	my $bevbu = $filename->get('burstgti', 'bat', '', 0);
    if(! -e $bevbu ){
       $log->error(1, "$bevbu file 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";

	my ($t,$err,$r,$max1,$max2,$min1,$min2,$maxerr,$minerr)=0;
	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  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;
	}
} #end of sub gif_light_curve