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.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.12 $
#
#
##############################################################################


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', 'evsh{ps,sl,as}', '*');
	unless(@eventFiles) {
		$log->entry("No filtered BAT event mode data\n");
		return;
	}

	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);

		$log->entry("TOTDUR=$info{totdur} T90=$info{t90dur} T50=$info{t50dur}");

		# 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);
	{
		$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", {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();

	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('T50DUR', $t50dur)  # T50 value and error
			->readkey('T50ERR', $t50err)
			->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;