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.27 2006/05/06 21:16:37 apsop
# HISTORY: Set enable_map to none if one is not found.
# HISTORY:
# 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.27 $
#
#
##############################################################################
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;
}
$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");
$enable_map = 'none';
}
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;