package Subs::BATIntervals;
##############################################################################
#
# DESCRIPTION: Compute the preliminary 4 ms light curve and the Bayesian
# blocks GTIS. This also gives the duration measures.
#
# Port of Craig Markwardt's BAT::bayes module to SDC.
#
# HISTORY: $Log: BATIntervals.pm,v $
# HISTORY: Revision 1.17 2007/10/12 20:32:31 apsop
# HISTORY: When running battblocks, do not treat warnings as errors. Do not read T50 values from output file, as they are sometimes missing and never used.
# HISTORY:
# HISTORY: Revision 1.16 2007/10/05 14:16:23 apsop
# HISTORY: Use all of the bat event data in the light curve for determining bat timing intervals.
# HISTORY:
# HISTORY: Revision 1.15 2007/09/12 18:04:58 apsop
# HISTORY: Fix bug in log entry call.
# HISTORY:
# HISTORY: Revision 1.14 2007/09/11 18:20:33 apsop
# HISTORY: Adjust trigger time for a short burst based on the FOREEXPO keyword.
# HISTORY:
# HISTORY: Revision 1.13 2007/01/31 21:16:41 apsop
# HISTORY: Do not do interval analysis on slew-only Bat data.
# HISTORY:
# HISTORY: Revision 1.12 2006/04/27 16:24:07 apsop
# HISTORY: Set seriousness for battblocks task to level 1, as the failure of this task is not necessarily an indication of a problem.
# HISTORY:
# HISTORY: Revision 1.11 2006/03/13 17:44:43 apsop
# HISTORY: Check for error in battblocks before proceeding. Reduces number of error messages when trying different binnings.
# HISTORY:
# HISTORY: Revision 1.10 2005/11/16 01:05:09 apsop
# HISTORY: Clean up intermediate gti files.
# HISTORY:
# HISTORY: Revision 1.9 2005/11/08 17:22:59 apsop
# HISTORY: Change calls for qualcal filenames. Clean up leftover temp files.
# HISTORY:
# HISTORY: Revision 1.8 2005/07/06 20:36:55 apsop
# HISTORY: Implemented BAT position refinement.
# HISTORY:
# HISTORY: Revision 1.7 2005/07/05 14:56:17 apsop
# HISTORY: Caught SDC up to date with BAT pipeline with exception of BAT
# HISTORY: position refinement.
# HISTORY:
# HISTORY: Revision 1.6 2005/06/01 16:03:07 apsop
# HISTORY: Change the burst gti file type from gti to burstgti.Subs/BATIntervals.pm
# HISTORY:
# HISTORY: Revision 1.5 2005/02/09 23:39:01 apsop
# HISTORY: Added message identifiers.
# HISTORY:
# HISTORY: Revision 1.4 2004/11/16 14:24:26 apsop
# HISTORY: Only try to process files that exist.
# HISTORY:
# HISTORY: Revision 1.3 2004/11/09 22:22:29 apsop
# HISTORY: Corrected name of peak GTI.
# HISTORY:
# HISTORY: Revision 1.2 2004/11/09 22:00:42 apsop
# HISTORY: Collect BAT GTIs in a single file.
# HISTORY:
# HISTORY: Revision 1.1 2004/11/08 18:40:09 apsop
# HISTORY: Module to calculate intermediate BAT GTIs.
# HISTORY:
# HISTORY: Revision 1.1 2004/09/24 20:52:21 wiegand
# HISTORY: Initial revision
# HISTORY:
#
# VERSION: $Revision: 1.17 $
#
#
##############################################################################
use Subs::Sub;
@ISA = ("Subs::Sub");
use strict;
use Util::HEAdas;
use Util::BATCave;
use Util::SwiftTags;
use Util::SimpleFITS;
use Astro::FITS::CFITSIO qw(:constants);
sub new {
my $proto=shift;
my $self=$proto->SUPER::new();
$self->{DESCRIPTION} = 'BAT interesting time interval determination';
return $self;
}
##################
# METHODS:
##################
sub body {
my $self = shift;
my $log = $self->log;
my $filename = $self->filename;
my $procpar = $self->procpar;
my $jobpar = $self->jobpar;
###############################################
# collect the event files
###############################################
my @eventFiles = $filename->getExisting('unfiltered', 'bat', 'evshsp');
unless(@eventFiles) {
$log->entry("No unfiltered BAT event mode data\n");
return;
}
##########################################################
# Bail if this is slew only data not from a burst
##########################################################
my $seq = $jobpar->read('sequence');
return
if( $seq%1000 != 0 && $seq%1000 < 899 && @eventFiles == 1 && $eventFiles[0] =~ /evshsl/ );
my %info = (
bblc => 'bb.lc',
inlist => join(',', @eventFiles),
qmap => $filename->get('qualcal', 'bat', ''),
bbfile => $filename->get('gti', 'bat', 'bayes', 0),
durfile => $filename->get('gti', 'bat', 'dur', 0),
peakint => 1,
tlookback => 10,
);
my $goodt90dur = 0;
my $goodtbinsize = 0;
foreach my $tbinsize (0.004, 0.064, 1.0, 16.0) {
$log->entry("trying bin size $tbinsize [s]");
$info{tbinsize} = $tbinsize;
$self->bayesCompute(\%info);
my $accepted = 'no';
if (not $info{complete}) {
# ignore results
}
else {
$info{complete} = 0;
my $t90dur = $info{t90dur};
if ($t90dur
and ($t90dur > 2.1 * $tbinsize)
and ($t90dur > 1.3 * $goodt90dur)) {
$goodt90dur = $t90dur;
$goodtbinsize = $tbinsize;
$accepted = 'yes';
}
$log->entry("totdur=$info{totdur} t90=$t90dur accepted=$accepted");
}
}
unlink($info{bblc});
my @toCopy; # files holding GTIs to copy
if ($goodt90dur == 0) {
# take from TDRSS message
unlink($info{bbfile}, $info{durfile});
my $tdrss_ack = Util::BATCave::get_tdrss_ack($self);
my $tdrss_gti = $filename->get('gti', 'bat', 'tdrss', 0);
if (-f $tdrss_ack) {
$log->error([ 1, BAT_TDRSS_GTIS ],
"BAT Bayesian block computation failed, will use TDRSS intervals");
$self->writeTDRSS_GTIs($tdrss_ack, $tdrss_gti);
@toCopy = ($tdrss_gti);
}
else {
$log->error([ 2, BAT_NO_GTIS ],
"Bayesian block computation failed, _and_ no TDRSS message");
}
}
else {
$log->entry("burst interval found (time bin size = $goodtbinsize)");
$log->entry("recomputing Bayesian blocks with best binning");
$info{tbinsize} = $goodtbinsize;
$self->bayesCompute(\%info);
my $durs = "TOTDUR=$info{totdur} T90=$info{t90dur}";
$durs .= " T50=$info{t50dur}" if $info{t50dur};
$log->entry($durs);
# Craig's code to read in the T90 and T50 values for the
# burst advocate report not implemented.
@toCopy = ($info{bbfile}, $info{durfile});
}
unlink($info{bblc});
# copy GTIs to the general file
my $gti_file = $filename->get('burstgti', 'bat', '', 0);
if (not $gti_file) {
$log->error([ 2, BAD_GTI ], "huh?! where did the BAT GTI file go?");
}
my $appender = Util::HEAdas->new('ftappend');
foreach my $file (@toCopy) {
my $fits = Util::SimpleFITS->open("<$file");
my $nhdu = $fits->nhdu;
$fits->close;
foreach my $ext (1 .. $nhdu-1) {
$appender->params({
infile => $file . "[$ext]",
outfile => $gti_file,
})
->run;
}
}
unlink @toCopy;
} # end of body method
sub writeTDRSS_GTIs
{
my ($self, $tdrss_ack, $tdrss_gti) = @_;
my $log = $self->log;
my $filename = $self->filename;
my $procpar = $self->procpar;
my $jobpar = $self->jobpar;
###############################################
# write fallback TDRSS intervals
###############################################
my ($mjdrefi, $mjdreff, $timezero);
my ($trigtime, $trigstop, $backstop, $trigsatf, $foreexpo);
{
$timezero = 0;
my $status = Util::SimpleFITS->open("<$tdrss_ack")
->readkey(MJDREFI => $mjdrefi)
->readkey(MJDREFF => $mjdreff)
# ->readkey(TIMEZERO => $timezero) # TODO: TIMEZERO is missing
->readkey(TRIGTIME => $trigtime)
->readkey(TRIGSTOP => $trigstop)
->readkey(BACKSTOP => $backstop)
->readkey(TRIGSATF => $trigsatf)
->readkey(FOREEXPO => $foreexpo)
->close
->status;
if ($status) {
$log->error([ 1, SIMPLEFITS_ERROR ],
"unable to read keys from BAT TDRSS alert message [$tdrss_ack]");
}
}
if($trigsatf >= 10000 && $trigsatf <= 10999){
$log->entry("Adjust trigger start and stop times for a short burst.");
$trigstop = $trigtime + 0.320 + $foreexpo/1000.0;
$trigtime = $trigtime - $foreexpo/1000.0;
}
my @gtis = (
{ name => 'GTI_TOT',
start => $trigtime,
stop => $trigstop,
},
{ name => 'GTI_BKG1',
start => -1e307,
stop => $backstop,
},
{ name => 'GTI_BKG2',
start => $trigstop + 60,
stop => 1e307,
},
);
# when using TDRSS intervals, fall back to GTI_TOT
foreach my $ext (qw(GTI_T90 GTI_T50 GTI_PEAK)) {
push(@gtis, {
name => $ext,
start => $trigtime,
stop => $trigstop,
comment => 'Bayesian block computation failed, using TDRSS intervals',
});
}
my $fits = Util::SimpleFITS->open("+<$tdrss_gti");
if (not $fits or $fits->status) {
$log->error([ 1, SIMPLEFITS_ERROR ],
"unable to create BAT TDRSS GTI file $tdrss_gti");
}
foreach my $gti (@gtis) {
my $status = $fits
->createtab($gti->{name})
->writekey(NAXIS2 => 1)
->insertcol({ TTYPE => [ "START", "GTI Start Time"],
TFORM => "1D", TUNIT => "s"})
->insertcol({ TTYPE => [ "STOP", "GTI Stop Time"],
TFORM => "1D", TUNIT => "s"})
->writecol("START", {type => TDOUBLE}, $gti->{start})
->writecol("STOP", {type => TDOUBLE}, $gti->{stop})
->status;
if ($status) {
$log->error([ 1, SIMPLEFITS_ERROR ],
"unable to create BAT TDRSS $gti->{name}");
}
if ($gti->{comment}) {
$fits->handle->write_comment($gti->{comment}, $status);
}
}
$fits->close;
undef($fits);
}
sub bayesCompute
{
my ($self, $info) = @_;
my $log = $self->log;
my $filename = $self->filename;
my $procpar = $self->procpar;
my $jobpar = $self->jobpar;
unlink($info->{bblc}, $info->{bbfile}, $info->{durfile});
my $batbinevt = Util::HEAdas->new('batbinevt')
->params({
infile => $info->{inlist},
outfile => $info->{bblc},
outtype => 'LC',
timebinalg => 'uniform',
energybins => '15-350', # keV # TODO: make parameter?
timedel => $info->{tbinsize},
detmask => $info->{qmap},
ecol => 'ENERGY',
weighted => 'YES',
outunits => 'RATE',
clobber => 'YES',
chatter => 2,
})
->run;
return if $batbinevt->had_error;
my $battblocks = Util::HEAdas->new('battblocks')
->params({
infile => $info->{bblc},
outfile => $info->{bbfile},
durfile => $info->{durfile},
tlookback => $info->{tlookback},
tpeak => $info->{peakint},
bkgsub => 'NO',
timecol => 'TIME',
countscol => 'RATE',
errcol => 'ERROR',
hduclas3 => 'RATE',
clobber => 'yes',
chatter => 2,
});
$battblocks->seriousness(1);
$battblocks->run();
return if $battblocks->had_error()==2;
if($battblocks->had_error() == 1){
my $errout = $battblocks->stderr();
unless( $errout =~ /^WARNING/m ){
return;
}
}
if (not -f $info->{bbfile}) {
$log->entry("no GTI found for timedel $info->{tbinsize}");
return;
}
# Read the Bayesian blocks...
# ... and also read the T50/90 keywords
my ($t90dur, $t90err, $t50dur, $t50err, $tpeak, $cntflu);
my ($totdur, $totstart, $totstop);
my @bbstart;
my @bbstop;
my $status = 0;
$status = Util::SimpleFITS->open("<$info->{bbfile}")
->move('STDGTI')
->readcol('START', TDOUBLE, [], \@bbstart)
->readcol('STOP', TDOUBLE, [], \@bbstop)
->readkey('T90DUR', $t90dur) # T90 value and error
->readkey('T90ERR', $t90err)
->readkey('TPEAK', $tpeak) # Epoch of peak flux
->readkey('CNTFLU', $cntflu) # Counts fluence for GRB
->readkey('TOTDUR', $totdur) # Total duration
->readkey('TOTSTART', $totstart) # Start of total burst
->readkey('TOTSTOP', $totstop)# End of total burst
->close()
->status();
if ($status) {
$log->error([ 1, SIMPLEFITS_ERROR ],
"could not open/read $info->{bbfile} [$status]");
return;
}
if (@bbstart < 3) {
$log->entry("warning: too few Bayesian blocks to bracket the GRB");
return;
}
$info->{t90dur} = $t90dur;
$info->{t50dur} = $t50dur;
$info->{totdur} = $totdur;
# Same as TOTSTART and TOTSTOP (in build 9.1)
my $tstart = $bbstop[0];
my $tstop = $bbstart[-1];
my $tmpfile = 'gti.tmp';
unlink($tmpfile);
my $copier = Util::HEAdas->new('ftcopy');
$copier->params({
infile => $info->{durfile} . '[1]',
outfile => $tmpfile,
copyall => 'no',
})
->run;
if ($copier->had_error or not -f $tmpfile) {
$log->error([ 2, HEATOOL_ERROR ], "could not create GTI_BKG1");
return;
}
# Write first background interval
$status = Util::SimpleFITS->open("+<$tmpfile")
->move(2)
->writecol("START", {type => TDOUBLE}, -1.0e307)
->writecol("STOP", {type => TDOUBLE}, $tstart)
->writekey("TSTART", -1.0e307)
->writekey("TSTOP", $tstart)
->writekey("EXTNAME", "GTI_BKG1")
->close()
->status();
if ($status) {
$log->error([ 2, HEATOOL_ERROR ],
"could not write GTI_BKG1 (status=$status)");
return;
}
my $appender = Util::HEAdas->new('ftappend');
$appender->params({
infile => $tmpfile . '[GTI_BKG1]',
outfile => $info->{durfile},
})
->run;
return if $appender->had_error;
# Write second background interval
$status = Util::SimpleFITS->open("+<$tmpfile")
->move(2)
->writecol("START", {type => TDOUBLE}, $tstop)
->writecol("STOP", {type => TDOUBLE}, 1.0e307)
->writekey("TSTART", $tstop)
->writekey("TSTOP", 1.0e307)
->writekey("EXTNAME", "GTI_BKG2")
->close()
->status();
if ($status) {
$log->error([ 2, SIMPLEFITS_ERROR ],
"could not write GTI_BKG2 (status=$status)");
return;
}
$appender->params({
infile => $tmpfile . '[GTI_BKG2]',
outfile => $info->{durfile},
})
->run;
return if $appender->had_error;
# Check for GTI_TOT (build 9.1 has it already)
my $fits = Util::SimpleFITS->open("<$info->{durfile}")
->move('GTI_TOT');
$status = $fits->status();
$fits->setstatus(0)->close();
undef($fits);
# If it does not exist, then create one
if ($status) {
$status = Util::SimpleFITS->open("+<$tmpfile")
->move(2)
->writecol("START",{type => TDOUBLE},$tstart)
->writecol("STOP",{type => TDOUBLE},$tstop)
->writekey("TSTART",$tstart)
->writekey("TSTOP",$tstop)
->writekey("EXTNAME","GTI_TOT")
->close()
->status();
if ($status) {
$log->error([ 2, SIMPLEFITS_ERROR ],
"could not write GTI_TOT [status=$status]");
return;
}
$appender->params({
infile => $tmpfile . '[GTI_TOT]',
outfile => $info->{durfile},
})
->run;
return if $appender->had_error;
}
unlink($tmpfile);
$info->{complete} = 1;
}
1;