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.26 2016/08/02 17:11:45 apsop
# HISTORY: Fix calculation of BAT GTIs for low significance triggers using scaled map
# HISTORY:
# HISTORY: Revision 1.25 2016/07/21 14:20:09 apsop
# HISTORY: Updates for 3.17.05. CALDB patches for XRT 20160609 and clock v114. xrtwtcorr v0.2.2. Fix issue retrieving BAT scaled map.
# HISTORY:
# HISTORY: Revision 1.24 2013/08/19 20:29:24 apsop
# HISTORY: Only log fetch_from_repository result if fetch succeeded, to
# HISTORY: avoid a Perl message about using an undefined variable. (If
# HISTORY: unsuccessful, fetch_from_repository already logs an error.)
# HISTORY:
# HISTORY: Revision 1.23 2013/05/30 07:54:54 apsop
# HISTORY: Log fetches from BAT repository
# HISTORY:
# HISTORY: Revision 1.22 2009/07/16 19:36:57 apsop
# HISTORY: in battblocks changed tstart, tstop to its real values
# HISTORY:
# HISTORY: Revision 1.21 2009/03/05 21:24:00 apsop
# HISTORY: change start=0 stop=1E100
# HISTORY:
# HISTORY: Revision 1.20 2008/12/10 13:52:22 apsop
# HISTORY: Fix problem with multiple copies of GTI extensions.
# HISTORY:
# HISTORY: Revision 1.19 2008/05/16 13:43:42 apsop
# HISTORY: Implement best algorithm for short burst trigger interval.
# HISTORY:
# HISTORY: Revision 1.18 2007/12/18 19:54:41 apsop
# HISTORY: Set lookback to 0 in battblock if bin size > 0.5 .
# HISTORY:
# 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.26 $
#
#
##############################################################################
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/ );
###################################
# Return if no burst alert message
###################################
my $tdrss_ack = Util::BATCave::get_tdrss_ack($self);
return unless ($tdrss_ack && -f $tdrss_ack);
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,
);
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;
$info{tlookback} = $tbinsize >= 0.5 ? 0 : 10;
$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);
$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->entry("burst interval found (time bin size = $goodtbinsize)");
$log->entry("recomputing Bayesian blocks with best binning");
$info{tbinsize} = $goodtbinsize;
$info{tlookback} = $goodtbinsize >= 0.5 ? 0 : 10;
$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.");
#########################
# Check for a scaled map
#########################
my $seq = $jobpar->read('sequence');
my $index = Util::BATCave::get_bscalemap_index($log, $seq);
$log->entry("BATIntervals.pm requesting from repository " .
"bscalemap: $index");
my $bscale = $filename->fetch_from_repository('bscalemap',
'bat', '', $index);
my $usedBscale = undef;
BSCALE_BAIL: {
if ($bscale && -f $bscale) {
$log->entry("fetched from repository " .
"bscalemap: \$bscale=$bscale");
$log->entry("Use $bscale to determine trigger times.");
=pod
-----Original Message-----
From: Markwardt, Craig B (GSFC-6620)
Sent: Friday, July 22, 2016 10:35 AM
Subject: Re: BAT shared repository bug
[snip]
I think the only information we have left to select on is the trigger score.
Bob, I think the procedure should be,
1. Open scaled map file.
2. Select only the rows where TIME is within +/- 1 sec of the known trigger time
3. Of those, select the row with maximum bestTriggerScore.
=cut
my $bscfits = Util::FITSfile->new($bscale, 1);
my @cols = qw(TIME
mStartForeSec mStartForeSubsec
mEndForeSec mEndForeSubsec
mEndBackSec mEndBackSubsec
bestTriggerScore);
# hack: Util::FITSfile limits length of columns in a broken way
$bscfits->{MAX_COL_CHARACTERS} = 10 + length(join(' ', @cols));
$bscfits->cols(@cols);
my @values = $bscfits->table();
if (@values % @cols) {
my $v = @values;
my $c = @cols;
$log->error(2, "loading bscalemap table failed [$v values for $c cols]");
last BSCALE_BAIL;
}
my @entries;
while (@values) {
# create a record for the next row
my %entry = map { ($_ => shift(@values)) } @cols;
push(@entries, \%entry);
}
# [2] and [3] are done in this loop
my $best = undef;
foreach my $e (@entries) {
if (abs($e->{TIME} - $trigtime) <= 1) {
if (not $best or $e->{bestTriggerScore} > $best->{bestTriggerScore}) {
$best = $e;
$log->entry("best so far is $e->{TIME} with score $e->{bestTriggerScore}");
}
}
}
if ($best) {
# m*Subsec has already been scaled to seconds (by 2e-5)
$trigtime = $best->{mStartForeSec} + $best->{mStartForeSubsec};
$trigstop = $best->{mEndForeSec} + $best->{mEndForeSubsec};
$backstop = $best->{mEndBackSec} + $best->{mEndBackSubsec};
$usedBscale = 1;
}
}
} # end BSCALE_BAIL
if (not $usedBscale) {
$log->entry("using defaults for trigtime and trigstop");
# foreexpo has already been scaled to seconds
$trigstop = $trigtime + 0.320 + $foreexpo;
$trigtime = $trigtime - $foreexpo;
}
}
# Check for partially sane values - an image trigger will
# be missing the BACKSTRT/STOP values
if ($trigtime != 0.0 && $trigstop != 0.0 && $backstop == 0.0) {
$backstop = $trigtime;
}
my @gtis = (
{ name => 'GTI_TOT',
start => $trigtime,
stop => $trigstop,
},
{ name => 'GTI_BKG1',
start => 0,
stop => $backstop,
},
{ name => 'GTI_BKG2',
start => $trigstop + 60,
stop => 1e100
},
);
# 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 $file=$info->{bblc};
my $fits=Util::FITSfile->new($file);
my $tstart=$fits->keyword('TSTART');
my $tstop =$fits->keyword('TSTOP');
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,
global_tstart => $tstart,
global_tstop => $tstop
});
$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;
$info->{complete} = 1;
}
1;