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