package Subs::BATBurst;
##############################################################################
#
# DESCRIPTION: This subroutine does calibration of BAT event data.
# DESCRIPTION: It runs bateconvert to do energy scale calibrations,
# DESCRIPTION: and then runs batmaskwtevt to calculate mask weights for
# DESCRIPTION: each event which correspond to the burst position.
#
# DESCRIPTION: Port of Craig Markwardt's BAT::maskwt module to SDC.
#
#
# HISTORY:
# HISTORY: $Log: BATBurst.pm,v $
# HISTORY: Revision 1.26 2005/03/07 20:22:31 apsop
# HISTORY: Comment out filtering of TIME<1000. This now done in StartEndTimes.
# HISTORY:
# HISTORY: Revision 1.25 2005/02/09 23:39:01 apsop
# HISTORY: Added message identifiers.
# HISTORY:
# HISTORY: Revision 1.24 2005/02/08 18:26:52 apsop
# HISTORY: Fix bateconvert call to work with build 11. Remove events with time of less than 1000.
# HISTORY:
# HISTORY: Revision 1.23 2004/12/07 17:34:57 apsop
# HISTORY: Change log message to error message when bat files cannot be split. Downstream error messages that result have been changed to log messages.
# HISTORY:
# HISTORY: Revision 1.22 2004/12/05 23:26:54 apsop
# HISTORY: Split event list after calibration
# HISTORY:
# HISTORY: Revision 1.21 2004/11/19 21:46:48 apsop
# HISTORY: New version of xrt2fits.
# HISTORY:
# HISTORY: Revision 1.20 2004/11/10 18:01:06 apsop
# HISTORY: Replace SimpleFITS call for getting TSTART with FITSfile call.
# HISTORY:
# HISTORY: Revision 1.19 2004/11/09 22:03:14 apsop
# HISTORY: Reverted mode of short events GTI.
# HISTORY:
# HISTORY: Revision 1.18 2004/11/08 18:38:40 apsop
# HISTORY: Updated with Craig's BAT::maskwt code.
# HISTORY:
# HISTORY: Revision 1.17 2004/10/25 14:47:27 apsop
# HISTORY: Fix bugs in splitting of DPH files.
# HISTORY:
# HISTORY: Revision 1.16 2004/10/24 20:00:55 apsop
# HISTORY: Rework of split_events method. Now make pemanent GTI file.
# HISTORY:
# HISTORY: Revision 1.15 2004/08/30 13:23:06 apsop
# HISTORY: Fix up parameters for batmaskwtevt, including using bat_z and origin_z
# HISTORY:
# HISTORY: Revision 1.14 2004/08/16 15:23:09 apsop
# HISTORY: Turn on writing of HISTORY keywords.
# HISTORY:
# HISTORY: Revision 1.13 2004/08/13 14:28:31 apsop
# HISTORY: Test for case where slew times cannot be determined.
# HISTORY:
# HISTORY: Revision 1.12 2004/06/29 14:21:13 apsop
# HISTORY: Add residfile and pulserfile to bateconvert call.
# HISTORY:
# HISTORY: Revision 1.11 2004/06/11 19:49:48 apsop
# HISTORY: Handle case where there is no pointing time.
# HISTORY:
# HISTORY: Revision 1.10 2004/05/06 20:02:33 dah
# HISTORY: Add version number back into the header comments.
# HISTORY:
# HISTORY: Revision 1.9 2004/04/16 20:21:18 dah
# HISTORY: Begin using embedded history records
# HISTORY:
#
# VERSION: 0.0
#
#
##############################################################################
use Subs::Sub;
@ISA = ("Subs::Sub");
use strict;
use Util::SwiftTags;
sub new {
my $proto=shift;
my $self=$proto->SUPER::new();
$self->{DESCRIPTION}="BAT Burst Processing";
return $self;
}
##################
# METHODS:
##################
sub body {
my $self=shift;
my $log =$self->log();
my $filename=$self->filename();
my $procpar =$self->procpar();
my $jobpar =$self->jobpar();
################################################
# Events file times of less than 1000 are bogus
################################################
my @event_files = $filename->get('unfiltered', 'bat', 'ev??to', '*');
# my $copy = Util::HEAdas->new('ftcopy');
# my $filter = '[EVENTS][TIME > 1000]';
foreach my $unf (@event_files){
my ($inst, $mode, $index) = $filename->parse($unf, 'unfiltered');
my $new_name = $filename->get('unfiltered', 'bat', substr($mode,0,4).'sp', $index);
# $copy->params({infile => $unf . $filter,
# outfile => $new_name})
# ->run();
# unlink $unf;
rename $unf, $new_name;
}
$self->calibrate_events;
$self->mask_weight_events;
###################################################
# Split up the bat event list based on time period
###################################################
$self->split_events();
} # end of body method
################################################################
# Split up the bat event lists and dph's based on time period
###############################################################
sub split_events{
my $self=shift;
my $log =$self->log();
my $filename=$self->filename();
my $procpar =$self->procpar();
my $jobpar =$self->jobpar();
my $root = $filename->sequence_specific();
my @dph_files = $filename->get('rawdph', 'bat', 'svto', '*');
my @ack_files = $filename->get('tdmess', 'bat', 'ce', '*');
@ack_files = grep /^${root}/, @ack_files;
my @event_files = $filename->get('unfiltered', 'bat', 'evshsp', '*');
unless( @ack_files ){
$log->error([ 1, BAT_NO_CENTROID ],
'No burst ack message, so will not split BAT files.');
foreach my $raw (@dph_files){
my ($inst, $mode, $index) = $filename->parse($raw, 'dph');
my $new_name = $filename->get('rawdph', 'bat', 'svpb'.substr($mode,4), $index);
rename $raw, $new_name;
}
return;
}
####################################
# Do the splitting of events lists
#####################################
my $batack = Util::FITSfile->new($ack_files[0], 0);
$log->entry("Split BAT event lists based on end of slew time.");
# TODO: want to rename this from po => as, but po set in Attitude.pm
my $spo_gti = Util::FITSfile->new($filename->get('gti', 's', 'po', 0), 'STDGTI');
my @start = $spo_gti->cols('START')->table();
my @stop = $spo_gti->cols('STOP')->table();
my $trigtime = $batack->keyword('TRIGTIME');
#######################################################################
# Slewing senarios which are handled:
#
# point ----------- -------------
# | | | |
# slew ----- ---------- ---------
# *
# trigger
#
#
#
# point ----- -------------
# | | | |
# slew ----------- ---------- ---------
# *
# trigger
#
#
#
# point ------------------------------------------------
# |
# slew -----
# *
# trigger
#
###########################################################################
my ($slew_end, $slew_start, $preslew_start, $postslew_stop);
if(@start){
$preslew_start = ( grep($_ < $trigtime, @start) )[-1];
$slew_start = ( grep($_ > $trigtime, @stop) )[0];
$slew_end = ( grep($_ > $trigtime, @start) )[0];
$postslew_stop = ( grep($_ > $trigtime, @stop) )[1];
if( $slew_start && $slew_end && $slew_end < $slew_start ){
$preslew_start = $slew_end;
$slew_end = ( grep($_ > $trigtime, @start) )[1];
$postslew_stop = ( grep($_ > $trigtime, @stop) )[2];
}
}
$preslew_start = 0 unless $preslew_start;
$slew_start = 1.0E12 unless $slew_start;
$slew_end = 1.1E12 unless $slew_end;
$postslew_stop = 1.2E12 unless $postslew_stop;
$log->entry("BAT trigger time base on $ack_files[0] is $trigtime .");
$log->entry("Slew start/stop based on ACS settled flags: $slew_start / $slew_end .");
$log->entry("Pre/Post slew start/stop: $preslew_start / $postslew_stop .");
##############################
# Set up to split event files
##############################
my $gti_pre = $filename->get('gti', 'bat', '', 0);
my $gti_post = 'bat_post_gti.tmp';
my $gti_slew = 'bat_slew_gti.tmp';
my $datafile = 'bat_data.tmp';
my $cdfile = 'bat_cols.tmp';
my $hdfile = 'bat_head.tmp';
open COLS, ">$cdfile";
open HEAD, ">$hdfile";
print COLS "START 1D s\n";
print COLS "STOP 1D s\n";
print HEAD "EXTNAME = 'GTI_PRE ' / name of this binary table extension\n";
print HEAD "HDUCLASS= 'OGIP ' / format conforms to OGIP/GSFC conventions\n";
print HEAD "HDUCLAS1= 'GTI ' / Extension contains Good Time Intervals\n";
print HEAD "HDUCLAS2= 'STANDARD' / Extension contains Good Time Intervals\n";
print HEAD "TIMESYS = 'TT ' / time measured from\n";
print HEAD "MJDREFI = 51910 / MJD reference day\n";
print HEAD "MJDREFF = 7.428703700000000E-04 / MJD reference (fraction of day)\n";
print HEAD "CLOCKAPP= F / default\n";
print HEAD "TIMEUNIT= 's ' / unit for time keywords\n";
close COLS;
close HEAD;
my $table = Util::Ftool->new('fcreate')
->params({cdfile => $cdfile,
datafile => $datafile,
headfile => $hdfile});
open DATA, ">$datafile";
print DATA "$preslew_start $slew_start\n";
close DATA;
$table->params({outfile => $gti_pre})
->run();
open DATA, ">$datafile";
print DATA "$slew_start $slew_end\n";
close DATA;
$table->params({outfile => $gti_slew})
->run();
open DATA, ">$datafile";
print DATA "$slew_end $postslew_stop\n";
close DATA;
$table->params({outfile => $gti_post})
->run();
my $extractor = Util::Extractor->new();
foreach my $unf (@event_files) {
$log->entry("Splitting file $unf.");
my ($inst, $mode, $index) = $filename->parse($unf, 'unfiltered');
my $submode = substr($mode, 0, 4);
$extractor->instrument('bat', 'EVENT')
->infiles($unf);
my $post_file = $filename->get('unfiltered', 'b', $submode.'as', $index);
$extractor->gti($gti_post)
->outfile($post_file, 'event')
->run();
unlink $post_file
unless Util::FITSfile->new($post_file, 'EVENTS')->nrows();
my $pre_file = $filename->get('unfiltered', 'b', $submode.'ps', $index);
$extractor->gti($gti_pre)
->outfile($pre_file, 'event')
->run();
unlink $pre_file
unless Util::FITSfile->new($pre_file, 'EVENTS')->nrows();
my $slew_file = $filename->get('unfiltered', 'b', $submode.'sl', $index);
$extractor->gti($gti_slew)
->outfile($slew_file, 'event')
->run();
unlink $slew_file
unless Util::FITSfile->new($slew_file, 'EVENTS')->nrows();
}
##############################################
# Combine gtis together into gti product file
##############################################
my $fappend = Util::Ftool->new('fappend')
->params({outfile => $gti_pre});
$fappend->params({infile => "$gti_slew\[1\]\[col EXTNAME='GTI_SLEW'\]"})
->run();
$fappend->params({infile => "$gti_post\[1\]\[col EXTNAME='GTI_POST'\]"})
->run();
unlink $gti_post, $gti_slew;
############################
# Set up to split DPH files
############################
my $gti_merged = 'gti_merged.tmp';
$fappend->params({infile => $gti_merged.'[1]'});
my $mgtime = Util::Ftool->new('mgtime')
->params({merge => 'AND'});
my $fdelhdu = Util::Ftool->new('fdelhdu')
->params({confirm => 'no',
proceed => 'yes'});
$log->entry('Do splitting of DPH files based on slew times.');
$gti_pre = 'bat_pre_gti.tmp';
open DATA, ">$datafile";
print DATA "0 $trigtime\n";
close DATA;
$table->params({outfile => $gti_pre})
->run();
open DATA, ">$datafile";
print DATA "$trigtime 1.E12\n";
close DATA;
$table->params({outfile => $gti_post})
->run();
unlink $datafile;
unlink $cdfile;
unlink $hdfile;
#################################################################
# Do the splitting of dphs
#################################################################
foreach my $unf (@dph_files) {
$log->entry("Splitting file $unf.");
my ($inst, $mode, $index) = $filename->parse($unf, 'dph');
my $gaoff = substr($mode, 4);
##########################
# Get the after burst part
##########################
my $post_file = $filename->get('rawdph', 'b', 'svab'.$gaoff, $index);
Util::Ftool->new('fcopy')
->params({infile => "$unf\[BAT_DPH\][gtifilter(\'${gti_post}\[1\]\')]",
outfile => $post_file})
->run();
my $post_fits = Util::FITSfile->new($post_file, 'BAT_DPH');
if( $post_fits->nrows() ){
$post_fits->keyword('TSTART', $trigtime);
$fdelhdu->params({infile => $post_file.'[GTI]'})
->run();
$mgtime->params({ingtis => "$unf\[GTI\], $gti_post\[1\]",
outgti => $gti_merged})
->run();
$fappend->params({outfile => $post_file})
->run();
$post_fits->ext('STDGTI');
$post_fits->keyword('TSTART', $trigtime);
$post_fits->keyword('EXTNAME', 'GTI');
unlink $gti_merged;
}else{
unlink $post_file;
}
########################
# Get the pre-burst part
########################
my $pre_file = $filename->get('rawdph', 'b', 'svpb'.$gaoff, $index);
Util::Ftool->new('fcopy')
->params({infile => "$unf\[BAT_DPH\][gtifilter(\'${gti_pre}\[1\]\')]",
outfile => $pre_file})
->run();
my $pre_fits = Util::FITSfile->new($pre_file, 'BAT_DPH');
if( $pre_fits->nrows() ){
$pre_fits->keyword('TSTOP', $trigtime);
$fdelhdu->params({infile => $pre_file.'[GTI]'})
->run();
$mgtime->params({ingtis => "$unf\[GTI\], $gti_pre\[1\]",
outgti => $gti_merged})
->run();
$fappend->params({outfile => $pre_file})
->run();
$pre_fits->ext('STDGTI');
$pre_fits->keyword('TSTOP', $trigtime);
$pre_fits->keyword('EXTNAME', 'GTI');
unlink $gti_merged;
}else{
unlink $pre_file;
}
unlink $unf;
}
unlink $gti_pre, $gti_post;
}
sub calibrate_events
{
my $self = shift;
my $log = $self->log;
my $filename = $self->filename;
my $procpar = $self->procpar;
my $jobpar = $self->jobpar;
###############################################
# make sure we have some event data to work on
###############################################
my @all_event_files = $filename->get('unfiltered', 'bat', 'evsh??', '*');
unless(@all_event_files) {
$log->entry("No BAT event mode data\n");
return;
}
############################################################
# fill calibrated energy column
# here we need to add a system to handle the calibration
# files. (shared across sequences)
###########################################################
$log->entry("Running bateconvert on all event files");
my %common = (
residfile => $filename->fetch_cal('resid', 'bat'),
pulserfile => $filename->fetch_cal('pulser', 'bat'),
fltpulserfile => $filename->fetch_cal('pulseflt', 'bat'),
outfile => 'NONE',
clobber => 'yes',
chatter => 3,
history => 'yes',
calmode => 'QUADRATIC',
zeroit => 'no',
scaled_energy => 'yes'
);
my $tool = Util::HEAdas->new('bateconvert')
->params(\%common);
foreach my $unf (@all_event_files) {
$log->entry("Calibrating $unf");
###############################################
# get the gain offset calibration file to use
# for this event file
###############################################
my $unf_fits = Util::FITSfile->new($unf, 0);
my $time = $unf_fits->keyword('TSTART');
unless($time) {
$log->error(2, "unable to calibrate $unf, missing TSTART");
next;
}
my $gainoff = $filename->fetch_from_repository('bgaoff', 'b', '', $time);
$log->entry("Using gain/offset cal file $gainoff based on a start time of $time .");
###########################
# now run the calibration
###########################
$tool->params({
calfile => $gainoff,
infile => $unf,
})
->run;
} # end of bateconvert loop over files
}
sub mask_weight_events
{
my $self = shift;
my $log = $self->log;
my $filename = $self->filename;
my $procpar = $self->procpar;
my $jobpar = $self->jobpar;
$log->entry('Computing BAT mask weighting for the desired source');
my @events = $filename->get('unfiltered', 'bat', 'evsh??', '*');
if (not @events) {
$log->entry('no BAT (short) event files to process');
return;
}
my %common = (
attitude => $filename->get('attitude', 's'),
ra => $jobpar->read('burst_ra'),
dec => $jobpar->read('burst_dec'),
bat_z => $jobpar->read('bat_z'),
aperture => $filename->fetch_cal('aperture', 'bat'),
detmask => $filename->get('qualcal', 'bat', 'cb'),
coord_type => 'sky',
rebalance => 'yes',
corrections => 'flatfield,ndets,pcode,maskwt',
teldef => $filename->fetch_cal('teldef', 'bat'),
origin_z => $jobpar->read('bat_origin_z'),
);
my $tool = Util::HEAdas->new('batmaskwtevt')
->params(\%common);
foreach my $file (@events) {
$tool->params({
infile => $file,
})
->run;
}
}
1;