Subs::BATSpectra (version $)


package Subs::BATSpectra;
##############################################################################
#
# DESCRIPTION: Accumulate BAT spectra and follow on products
#
#		Port of Craig Markwardt's BAT::spectra, BAT::lc and setup to SDC.
#
#
# HISTORY: $Log: BATSpectra.pm,v $
# HISTORY: Revision 1.26  2005/12/19 16:13:18  apsop
# HISTORY: Add date keywords to primary header of mask tagged lightcurves.
# HISTORY:
# HISTORY: Revision 1.25  2005/11/08 17:28:44  apsop
# HISTORY: Change calls for qualcal filenames. Comment out temporary fix to enable map extension name.
# HISTORY:
# HISTORY: Revision 1.24  2005/08/24 12:41:45  apsop
# HISTORY: Bugzilla SDC 274 [from Craig Markwardt].
# HISTORY: Summary: BAT GRB processing - make RATE spectra instead of COUNTS
# HISTORY: "The SDC pipeline has been defaulting to making COUNTS spectra
# HISTORY: for BAT GRB processing.
# HISTORY: "It turns out that the standard downstream processing tools like
# HISTORY: GRPPHA and CMPPHA are not very good at handling type II counts
# HISTORY: spectra.  As a workaround, I recommend that the SDC switch
# HISTORY: the BAT spectra to be RATE spectra."
# HISTORY:
# HISTORY: Revision 1.23  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.22  2005/06/01 16:02:19  apsop
# HISTORY: Change method for fetching quad rate files, as rate types have change.
# HISTORY:
# HISTORY: Revision 1.21  2005/05/16 14:15:33  apsop
# HISTORY: Handle missing input files more gracefully when creating images.  Pass
# HISTORY: CALDB for ebounds parameter when creating mask tagged light curves.
# HISTORY:
# HISTORY: Revision 1.20  2005/05/13 19:38:02  apsop
# HISTORY: Implemented product agreements between HEASARC and BAT team.
# HISTORY:
# HISTORY: Revision 1.19  2005/05/04 21:29:28  apsop
# HISTORY: Create mask tagged light curves.
# HISTORY:
# HISTORY: Revision 1.18  2005/04/08 01:15:14  apsop
# HISTORY: Remove unweighted 4ms lc
# HISTORY:
# HISTORY: Revision 1.17  2005/04/06 15:41:53  apsop
# HISTORY: Change to using CALDB for cal parameters.
# HISTORY:
# HISTORY: Revision 1.16  2004/12/31 01:36:26  apsop
# HISTORY: Comment out production of spectra which are not yet official.
# HISTORY:
# HISTORY: Revision 1.15  2004/12/06 01:00:12  apsop
# HISTORY: Test for presence of event files, change some error() calls than entry() calls.
# HISTORY:
# HISTORY: Revision 1.14  2004/11/29 22:14:35  apsop
# HISTORY: Corrected efile+escale parameters to batdrmgen.  Only process existing
# HISTORY: files.
# HISTORY:
# HISTORY: Revision 1.13  2004/11/16 14:30:25  apsop
# HISTORY: Avoid using undefined value.
# HISTORY:
# HISTORY: Revision 1.12  2004/11/09 22:22:29  apsop
# HISTORY: Corrected name of peak GTI.
# HISTORY:
# HISTORY: Revision 1.11  2004/11/09 22:00:42  apsop
# HISTORY: Collect BAT GTIs in a single file.
# HISTORY:
# HISTORY: Revision 1.10  2004/11/08 18:41:54  apsop
# HISTORY: Implemented Craig's Markwardt's BAT::spectra and BAT::lc modules.
# HISTORY:
#
# VERSION: $Revision: 1.26 $
#
#
##############################################################################


use Subs::Sub;

@ISA = ("Subs::Sub");
use strict;

use Util::HEAdas;
use Util::BATCave;


sub new {
	my $proto=shift;
	my $self=$proto->SUPER::new();

	$self->{DESCRIPTION}= 'Build BAT spectra, response matrices, light curves';

	return $self;
}

##################
# METHODS:
##################

sub body {

	my $self = shift;

	$self->make_spectra;

	$self->make_response_matrices;

	$self->make_light_curves;

	$self->make_mask_tagged_light_curves;

}



sub make_spectra {

	my $self = shift;

	my $log      = $self->log;
	my $filename = $self->filename;
	my $procpar  = $self->procpar;
	my $jobpar   = $self->jobpar;


	$log->entry('Accumulating spectra');

	my @events = $filename->getExisting('unfiltered', 'bat', 'evsh{ps,sl,as}');
	if (not @events) {
		$log->entry('no BAT preslew,slew,postslew events, not creating spectra');
		return;
	}

	# First the full-burst type intervals, use all the events
	my %info = (
		files => \@events,
		qmap  => $filename->get('qualcal', 'bat', ''),
	);


	###########################################################
	# only the preslew, slew and postslew spectra are HEASARC products
	###########################################################
	# $info{ebins} = Util::BATCave::chan4;
	# $info{ebins} = $filename->fetch_cal('ebounds', 'bat');
	$info{ebins} = 'CALDB';

	$info{outfile} = $filename->get('spectrum', 'bat', 'evps', 0);
	$info{gti} = Util::BATCave::get_gti($self, 'GTI_TOT');
	$self->makeSpectrum(\%info)
		if $info{gti};

	$info{outfile} = $filename->get('spectrum', 'bat', 'evsl', 0);
	$info{gti} = 'NONE';
	$self->makeSpectrum(\%info)
		if $info{gti};

	$info{outfile} = $filename->get('spectrum', 'bat', 'evas', 0);
	$info{gti} = 'NONE';
	$self->makeSpectrum(\%info)
		if $info{gti};

}


sub makeSpectrum
{
	my ($self, $info) = @_;

	my $log = $self->log;

	my $inlist = join(',', @{ $info->{files} });

	my $qmap = $info->{qmap} || 'NONE';
	my $gti = $info->{gti} || 'NONE';

	unlink($info->{outfile});

	##############################################
	# set up batbinevt
	#############################################
	my $binner = Util::HEAdas->new('batbinevt')
			->params({
					infile     => $inlist,
					outfile    => $info->{outfile},
					outtype    => 'PHA',
					timedel    => 0,
					timebinalg => 'uniform',
					energybins => $info->{ebins},
					detmask    => $qmap,
					ecol       => 'ENERGY',
					weighted   => 'YES',
					outunits   => 'RATE',
					gtifile    => $info->{gti},
					clobber    => 'YES',

					tstart     => 'INDEF',
					tstop      => 'INDEF',
					snrthresh  => 6.0,
					buffersize => 16384,
					chatter    => 5,
					history    => 'YES',
				})
			->run;

	####################
	# check for errors 
	####################
	if ($binner->had_error || ! -s $info->{outfile}) {
		$log->entry("Did not sucessfully make a spectrum for $inlist");
	}
}



sub make_response_matrices
{
	my $self = shift;

	my $log      = $self->log;
	my $filename = $self->filename;
	my $procpar  = $self->procpar;
	my $jobpar   = $self->jobpar;

	$log->entry('Computing response matrices');

	my $hkfile = $filename->get('bdaphk', 'bat') || 'NONE';

	my %common = (
		calfile => 'CALDB',
		depthfile => 'CALDB',

		escale => 'FILE',
		efile => 'CALDB',
		hkfile => $hkfile,

		chatter => 3,
		clobber => 'yes',
	);

	my $batdrmgen = Util::HEAdas->new('batdrmgen')
			->params(\%common);


	# preslew
	my $spectrum = $filename->get('spectrum', 'bat', 'evps');
	if (-f $spectrum) {
		$batdrmgen->params({
					infile => $spectrum,
					outfile => $filename->get('response', 'bat', 'evps', 0),
				})
			->run;
	}

	# postslew
	$spectrum = $filename->get('spectrum', 'bat', 'evas');
	if (-f $spectrum) {
		$batdrmgen->params({
					infile => $spectrum,
					outfile => $filename->get('response', 'bat', 'evas', 0),
				})
			->run;
	}


}



sub make_light_curves
{
	my $self = shift;

	my $log      = $self->log;
	my $filename = $self->filename;
	my $procpar  = $self->procpar;
	my $jobpar   = $self->jobpar;


	$log->entry('Accumulating light curves');

	my @events = $filename->getExisting('unfiltered', 'bat', 'evsh{ps,sl,as}');
	if (not @events) {
		$log->entry('no BAT preslew,slew,postslew events, not creating light curves');
		return;
	}

	my @config = (
		{ desc => 'weighted 1s light curve',
			outfile => $filename->get('lightcurve', 'bat', 'ev1s', 0),
			ebins => Util::BATCave::chan4,
			gti => 'NONE',
			tbinsize => 1.0,
			weighted => 'YES',
		},
		{ desc => 'weighted 64ms 1-channel light curve',
			outfile => $filename->get('lightcurve', 'bat', 'evms', 0),
			ebins => Util::BATCave::chan4,
			gti => 'NONE',
			tbinsize => 0.064,
			weighted => 'YES',
		},
	);

	my $qmap = $filename->get('qualcal', 'bat', '') || 'NONE';

	foreach my $c (@config) {
		$c->{files} = \@events;
		$c->{qmap} = $qmap;
		$self->lcMake($c);
	}

}


# BAT::lc->make

sub lcMake
{
	my ($self, $info) = @_;

	my $log      = $self->log;
	my $filename = $self->filename;
	my $procpar  = $self->procpar;
	my $jobpar   = $self->jobpar;


	my %par = (
		infile => join(',', @{ $info->{files} }),
		outfile => $info->{outfile},
		outtype => 'LC',

		energybins => $info->{ebins},
		detmask => $info->{qmap},
		weighted => $info->{weighted} || 'NO',
		outunits => 'RATE',
		ecol => 'ENERGY',
		gtifile => $info->{gti} || 'NONE',
		clobber => 'yes',
	);

	if ($info->{tbinsize} =~ /^g/i) {
		$par{timedel} = 0;
		$par{timebinalg} = 'gti';
	}
	else {
		$par{timedel} = $info->{tbinsize};
		$par{timebinalg} = 'uniform';
	}

	my $tool = Util::HEAdas->new('batbinevt')
			->params(\%par)
			->run;

}



sub make_mask_tagged_light_curves
{
	my $self = shift;

	my $log      = $self->log;
	my $filename = $self->filename;
	my $procpar  = $self->procpar;
	my $jobpar   = $self->jobpar;

	$log->entry('Making mask tagged light curves');

	my $enable_map = undef;

	my $tool = Util::HEAdas->new('batmasktaglc')
			->params({
				ebounds => 'CALDB',
			});

	foreach my $infile ($filename->get('rawmtlc', 'bat')) {

		my $start = 0;
		my $status = Util::SimpleFITS->open('<' . $infile)
			->readkey(TSTART => $start)
			->move('MASK_TAG_WEIGHT')
			->close
			->status;

		if ($status) {
			$log->entry("$infile does not have MASK_TAG_WEIGHTs");
			next;
		}

		# hack, apparently bat2fits writing Det_Enable_Map but
		# batmasktaglc expecting BAT_FLAGS
		$enable_map = $filename->fetch_from_repository('bdetflag', 'bat', '', $start);
#		if (not $enable_map) {
#			$log->error(1, "unable to get BAT enable map for time $start");
#			next;
#		}
#		rename $enable_map, 'enable_map.tmp';
#		$enable_map = 'enable_map.tmp';
#		$filename->fetch_from_repository('bdetflag', 'bat', '', $start);
#
#		Util::HEAdas->new('fthedit')
#		            ->params({
#				      infile => $enable_map . '[1]',
#				      keyword => 'EXTNAME',
#				      value => 'BAT_FLAGS',
#				      operation => 'add',
#				     })
#			     ->run;


		my $outfile = $infile;
		$outfile =~ s/_rw//;
		unlink($outfile)
			or $log->error(1, "unable to remove $outfile [$!]")
			if -e $outfile;

		my $quadfile = (grep /brtqd\.lc/, $filename->get('batrates', 'b', ''))[0];
		if (not $quadfile) {
			$log->error(1, "missing quadrate file for $infile");
			next;
		}

		$tool->params({
				infile    => $infile,
				quadfile  => $quadfile,
				maskwt    => $infile,
					# BAT2FITS appends the mask weight table to the raw LC
				outfile   => $outfile,
				detmask   => $enable_map,
			})
			->run;

		#######################################
		# Copy DATE keywords to primary header
		#######################################
		my $mtfits = Util::FITSfile->new($outfile);
		my $obs = $mtfits->keyword('DATE-OBS');
		my $end = $mtfits->keyword('DATE-END');
		$mtfits->ext(0);
		$mtfits->keyword('DATE-OBS', $obs);
		$mtfits->keyword('DATE-END', $end);
	}

}


1;