Subs::BATIntervals (version $)


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;