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


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;

}



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', 'cb'),
	);


#	$info{ebins} = $filename->fetch_cal('ebounds', 'bat');
#
#	$info{outfile} = $filename->get('spectrum', 'bat', 'evto', 0);
#	$info{gti} = Util::BATCave::get_gti($self, 'GTI_TOT');
#	$self->makeSpectrum(\%info)
#		if $info{gti};
#
#	$info{outfile} = $filename->get('spectrum', 'bat', 'evt9', 0);
#	$info{gti} = Util::BATCave::get_gti($self, 'GTI_T90');
#	$self->makeSpectrum(\%info)
#		if $info{gti};
#
#	$info{outfile} = $filename->get('spectrum', 'bat', 'evt5', 0);
#	$info{gti} = Util::BATCave::get_gti($self, 'GTI_T50');
#	$self->makeSpectrum(\%info)
#		if $info{gti};
#
#	$info{outfile} = $filename->get('spectrum', 'bat', 'evpk', 0);
#	$info{gti} = Util::BATCave::get_gti($self, 'GTI_PEAK');
#	$self->makeSpectrum(\%info)
#		if $info{gti};
#
#
	###########################################################
	# these are Craig's totflu and peakflu spectra...
	###########################################################
#	$info{ebins} = Util::BATCave::chan4;
#
#	$info{outfile} = $filename->get('spectrum', 'bat', 'evtf', 0);
#	$info{gti} = Util::BATCave::get_gti($self, 'GTI_TOT');
#	$self->makeSpectrum(\%info)
#		if $info{gti};
#
#	$info{outfile} = $filename->get('spectrum', 'bat', 'evpf', 0);
#	$info{gti} = Util::BATCave::get_gti($self, 'GTI_TOT');
#	$self->makeSpectrum(\%info)
#		if $info{gti};
#

	###########################################################
	# now the preslew, slew and postslew spectra
	###########################################################
	$info{ebins} = Util::BATCave::chan4;

	$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   => 'COUNTS',
					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 %common = (
		calfile => 'CALDB',
		depthfile => 'CALDB',

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

		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 4ms light curve',
			outfile => $filename->get('lightcurve', 'bat', 'evbu', 0),
			ebins => Util::BATCave::chan1,
			gti => 'NONE',
			tbinsize => 0.004,
			weighted => 'YES',
		},
		{ desc => 'weighted 64ms 1-channel light curve',
			outfile => $filename->get('lightcurve', 'bat', 'evms', 0),
			ebins => Util::BATCave::chan1,
			gti => 'NONE',
			tbinsize => 0.064,
			weighted => 'YES',
		},
#		{ desc => 'unweighted 4ms light curve',
#			outfile => $filename->get('rawlc', 'bat', 'evbu', 0),
#			ebins => Util::BATCave::chan1,
#			gti => 'NONE',
#			tbinsize => 0.004,
#			weighted => 'NO',
#		},
	);

	my $bbfile = Util::BATCave::get_gti($self, 'GTI_BAYES');

	push(@config,
		{ desc => 'weighted 4-channel light curve, Bayesian blocks',
			outfile => $filename->get('lightcurve', 'bat', 'evbb4c', 0),
			ebins => Util::BATCave::chan4,
			gti => $bbfile,
			tbinsize => 'g',
			weighted => 'YES',
		},
		{ desc => 'weighted 1-channel light curve, Bayesian blocks',
			outfile => $filename->get('lightcurve', 'bat', 'evbb1c', 0),
			ebins => Util::BATCave::chan1,
			gti => $bbfile,
			tbinsize => 'g',
			weighted => 'YES',
		},
		) if ($bbfile and -f $bbfile);

	my $qmap = $filename->get('qualcal', 'bat', 'cb') || '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;

}