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.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.5 $
#
#
##############################################################################
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
###############################################
# TODO: evshpo or as?
my @eventFiles = $filename->getExisting(
'unfiltered', 'bat', 'evsh{ps,sl,as}', '*');
unless(@eventFiles) {
$log->entry("No filtered BAT event mode data\n");
return;
}
my %info = (
lc4ms => 'lc4ms.tmp',
# lc4ms => $filename->get('???', 'bat', '*', 0), # TODO: keep it?
inlist => join(',', @eventFiles),
qmap => $filename->get('qualcal', 'bat', 'cb'),
bbfile => $filename->get('gti', 'bat', 'bayes', 0),
durfile => $filename->get('gti', 'bat', 'dur', 0),
);
foreach my $tbinsize (0.004, 0.064, 1.0, 16.0) {
$info{tbinsize} = $tbinsize;
$self->bayesCompute(\%info);
last if $info{complete};
}
unlink($info{lc4ms});
my @toCopy; # files holding GTIs to copy
if (not $info{complete}) {
# 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 {
@toCopy = ($info{bbfile}, $info{durfile});
}
# copy GTIs to the general file
my $gti_file = $filename->get('gti', 'bat');
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;
}
}
} # 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);
{
$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)
->close
->status;
if ($status) {
$log->error([ 1, SIMPLEFITS_ERROR ],
"unable to read keys from BAT TDRSS alert message [$tdrss_ack]");
}
}
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", [], $gti->{start}, TDOUBLE)
->writecol("STOP", [], $gti->{stop}, TDOUBLE)
->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->{lc4ms}, $info->{bbfile}, $info->{durfile});
my $batbinevt = Util::HEAdas->new('batbinevt')
->params({
infile => $info->{inlist},
outfile => $info->{lc4ms},
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->{lc4ms},
outfile => $info->{bbfile},
durfile => $info->{durfile},
tlookback => 10.0, # TODO: make parameter ?
tpeak => 1.0, # TODO: make parameter ?
bkgsub => 'NO',
timecol => 'TIME',
countscol => 'RATE',
errcol => 'ERROR',
hduclas3 => 'RATE',
clobber => 'yes',
chatter => 2,
})
->run;
# Read the Bayesian blocks...
# ... and also read the T50/90 keywords
my ($t90dur, $t90err, $t50dur, $t50err, $tpeak, $cntflu);
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('T50DUR', $t50dur) # T50 value and error
->readkey('T50ERR', $t50err)
->readkey('TPEAK', $tpeak) # Epoch of peak flux
->readkey('CNTFLU', $cntflu) # Counts fluence for GRB
# Build 9.1
# ->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;
}
# 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", [], -1.0e307, TDOUBLE)
->writecol("STOP", [], $tstart, TDOUBLE)
->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", [], $tstop, TDOUBLE)
->writecol("STOP", [], 1.0e307, TDOUBLE)
->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",[],$tstart,TDOUBLE)
->writecol("STOP",[],$tstop,TDOUBLE)
->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;