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.34 2007/02/05 15:27:58 apsop
# HISTORY: Trap errors from batmasktaglc that warn of no GTIs.
# HISTORY:
# HISTORY: Revision 1.33 2007/01/31 21:18:18 apsop
# HISTORY: Do not make products from slew-only event data.
# HISTORY:
# HISTORY: Revision 1.32 2006/10/26 21:27:31 apsop
# HISTORY: Apply GTI_TOT when extracting spectra. Trap no good time warning. Update spectra keywords.
# HISTORY:
# HISTORY: Revision 1.31 2006/09/28 18:31:59 apsop
# HISTORY: Fix bug in time intervals used for producing bat spectra, which were completely wrong.
# HISTORY:
# HISTORY: Revision 1.30 2006/07/25 19:49:34 apsop
# HISTORY: Add in calls to batupdatephakw and batphasyserr.
# HISTORY:
# HISTORY: Revision 1.29 2006/06/28 17:12:16 apsop
# HISTORY: Use new target pointing GTI for filtering mask weighted light curves.
# HISTORY:
# HISTORY: Revision 1.28 2006/06/15 22:21:36 apsop
# HISTORY: Remove "unsettled" time mask tagged lightcurves, and remove empty lightcurves.
# HISTORY:
# 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.34 $
#
#
##############################################################################
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;
my $log = $self->log;
my $filename = $self->filename;
my $procpar = $self->procpar;
my $jobpar = $self->jobpar;
##########################################################
# Bail if this is slew only data not from a burst
##########################################################
my @eventFiles = $filename->getExisting('unfiltered', 'bat', 'evsh{ps,sl,as}', '*');
my $seq = $jobpar->read('sequence');
return
if( $seq%1000 != 0 && $seq%1000 < 899 && @eventFiles == 1 && $eventFiles[0] =~ /evshsl/ );
$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;
my %info = (
ebins => 'CALDB',
gti => $filename->get('burstgti', 'bat', '', 0) . '[GTI_TOT]',
files => \@events,
qmap => $filename->get('qualcal', 'bat', ''),
auxfile => $filename->get('eventaux', 'b', '', 0)
);
###########################################################
# only the preslew, slew and postslew spectra are HEASARC products
###########################################################
@events = $filename->get('unfiltered', 'bat', 'evshps', '*');
if( @events ){
$info{outfile} = $filename->get('spectrum', 'bat', 'evps', 0);
$self->makeSpectrum(\%info)
}
@events = $filename->get('unfiltered', 'bat', 'evshas', '*');
if( @events ){
$info{outfile} = $filename->get('spectrum', 'bat', 'evas', 0);
$self->makeSpectrum(\%info)
}
@events = $filename->get('unfiltered', 'bat', 'evshsl', '*');
if( @events ){
$info{outfile} = $filename->get('spectrum', 'bat', 'evsl', 0);
$self->makeSpectrum(\%info)
}
}
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',
});
$binner->seriousness(1);
$binner->run();
if( $binner->had_error() == 2 ){
$log->error(2, "batbinevt level 2 error");
}elsif( $binner->had_error() == 1 ){
##################################################################
# Note: this case is not fatal, *IF* there were no overlapping
# good times between the input file and the good time interval
# file. In that case batbinevt does not create an output file.
#################################################################
my $errout = $binner->stderr();
foreach ( $errout =~ /^.*$/gm ){
if( ! /^=+$/ && ! /WARNING: no overlapping good time intervals were found/ ){
$log->error(2, "batbinevt level 2 error");
last;
}
}
}
####################
# check for errors
####################
if ($binner->had_error || ! -s $info->{outfile}) {
$log->entry("Did not sucessfully make a spectrum for $inlist");
return;
}
my $specfits = Util::FITSfile->new($info->{outfile}, 1);
my ($tstart, $tstop) = ( $specfits->keyword('TSTART'), $specfits->keyword('TSTOP') );
my $exp = $specfits->keyword('EXPOSURE');
my $start = Util::Date->new($tstart);
my $stop = Util::Date->new($tstop);
my $nhdus = $specfits->nhdus();
for(my $ext=0; $ext<$nhdus; $ext++){
next if $ext == 1;
$specfits->ext($ext);
$specfits->begin_many_keywords();
$specfits->keyword('TSTART', $tstart );
$specfits->keyword('TSTOP', $tstop );
$specfits->keyword('DATE-OBS', $start->date() .'T'.$start->time() );
$specfits->keyword('DATE-END', $stop->date() .'T'.$stop->time() );
$specfits->keyword('EXPOSURE', $exp );
$specfits->end_many_keywords();
}
if( -f $info->{auxfile} ){
Util::HEAdas->new('batupdatephakw')
->params({infile => $info->{outfile},
auxfile => $info->{auxfile}
})
->run();
}
Util::HEAdas->new('batphasyserr')
->params({infile => $info->{outfile},
syserrfile => 'CALDB'
})
->run();
}
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',
});
$tool->seriousness(1);
my $pointing = $filename->get('gti', 's', 'tp', 0);
my $tempwt = 'maskwt.tmp';
foreach my $infile ($filename->get('rawmtlc', 'bat')) {
###########################################################
# Filter out times when the s/c wasn't settled, and remove
# any resulting empty files
###########################################################
if( -f $pointing ){
my $copy = Util::HEAdas->new('ftcopy')
->params({infile => $infile . '[1][gtifilter(\"'. $pointing .'\")]',
outfile => $tempwt});
$copy->run();
my $newwt = Util::FITSfile->new($tempwt);
if($newwt->nrows() == 0){
$log->entry("Deleting $infile - empty after filtering.");
unlink $infile, $tempwt;
next;
}
my $stats=Util::HEAdas->new("ftstat")
->params({infile => $tempwt .'[1][col TIME]',
centroid => 'no'})
->verbose(0);
$stats->run();
my $parfile = $stats->parfile();
my ($pstart, $pstop) = ( $parfile->read('min'), $parfile->read('max') );
my $start = Util::Date->new( $pstart );
my $stop = Util::Date->new( $pstop );
foreach my $ext (0,1){
$newwt->ext($ext);
$newwt->begin_many_keywords();
$newwt->keyword('TSTART', $pstart );
$newwt->keyword('TSTOP', $pstop );
$newwt->keyword('DATE-OBS', $start->date() .'T'.$start->time() );
$newwt->keyword('DATE-END', $stop->date() .'T'.$stop->time() );
$newwt->end_many_keywords();
}
unlink $infile;
rename $tempwt, $infile;
}
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();
if( $tool->had_error() == 2 ){
$log->error(2, "batmasktaglc level 2 error");
}elsif( $tool->had_error() == 1 ){
my $errout = $tool->stderr();
foreach ( $errout =~ /^.*$/gm ){
unless( /WARNING: There were no GTIs. Is that right?/ ){
$log->error(2, "batmasktaglc level 2 error");
last;
}
}
}
#######################################
# 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;