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