Subs::BATBurst (version 0.0)


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.36  2006/04/27 16:27:02  apsop
# HISTORY: Set parameter pcodethresh to 0.05 in batmaskwtevt task.
# HISTORY:
# HISTORY: Revision 1.35  2006/01/16 17:50:21  apsop
# HISTORY: Dont split or update TSTART of dph files that dont need it.
# HISTORY:
# HISTORY: Revision 1.34  2005/12/19 16:10:32  apsop
# HISTORY: Add timing keywords to evaux file - should really be done in bat tool.
# HISTORY:
# HISTORY: Revision 1.33  2005/11/16 01:03:20  apsop
# HISTORY: Only make evaux file for slew-pointing event data.
# HISTORY:
# HISTORY: Revision 1.32  2005/11/08 17:01:20  apsop
# HISTORY: Change qualcal call for new filename.
# HISTORY:
# HISTORY: Revision 1.31  2005/07/05 14:56:16  apsop
# HISTORY: Caught SDC up to date with BAT pipeline with exception of BAT
# HISTORY: position refinement.
# HISTORY:
# HISTORY: Revision 1.30  2005/06/01 16:03:07  apsop
# HISTORY: Change the burst gti file type from gti to burstgti.Subs/BATIntervals.pm
# HISTORY:
# HISTORY: Revision 1.29  2005/04/19 15:40:18  apsop
# HISTORY: Change calmode parameter to INDEF.
# HISTORY:
# HISTORY: Revision 1.28  2005/04/06 15:40:31  apsop
# HISTORY: Change to using CALDB for cal parameters.
# HISTORY:
# HISTORY: Revision 1.27  2005/03/25 22:53:53  apsop
# HISTORY: Change pointing gti extname form STDGTI to GTI
# HISTORY:
# 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


sub determine_trigger_time
{
    my $self = shift;

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

    my $root = $filename->sequence_specific();

# this could be done once and stored in $jobpar

    my $ack = Util::BATCave::get_tdrss_ack($self);
 
    if (not $ack) {
      $log->error([ 1, BAT_NO_CENTROID ],
              'No burst ack message to provide trigger time.');
      return;
    }

    my $batack = Util::FITSfile->new($ack, 0);
    my $trigtime = $batack->keyword('TRIGTIME');

    $log->entry("BAT trigger time based on $ack is $trigtime .");
 
    return $trigtime;
}



################################################################
# 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 @dph_files = $filename->get('rawdph', 'bat', 'svto', '*');

    my @event_files = $filename->get('unfiltered', 'bat', 'evshsp', '*');

    my $trigtime = $self->determine_trigger_time;

    if (not $trigtime) {
      $log->error([ 1, BAT_NO_CENTROID ],
              'No trigger time so will not split BAT event 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
    #####################################

    $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 $gtipath =  $filename->get('gti', 's', 'po', 0);
    my $spo_gti;
    my @start;
    my @stop;

    if (-f $gtipath) {
       $spo_gti =  Util::FITSfile->new($gtipath, 'GTI');

       @start = $spo_gti->cols('START')->table();
       @stop =  $spo_gti->cols('STOP')->table();
    }
    else {
      $log->error([ 1, BAT_NO_GTIS ],
              'No spo GTI so will not split BAT event files.');
    }

    #######################################################################
    # 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("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('burstgti', '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);
      my $post_start = Util::FITSfile->new($unf)->keyword('TSTART');
      if($post_start >= $trigtime){
	rename $unf, $post_file;
	next;
      }

      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 => 'CALDB',
		pulserfile => 'CALDB',
		fltpulserfile => 'CALDB',
		outfile => 'NONE',
		clobber => 'yes',
		chatter => 3,
		history => 'yes',
		calmode => 'INDEF',
		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 => 'CALDB',
		detmask  => $filename->get('qualcal', 'bat', ''),
		coord_type => 'sky',
		rebalance => 'yes',
		corrections => 'flatfield,ndets,pcode,maskwt',
		teldef => 'CALDB',
		origin_z => $jobpar->read('bat_origin_z'),
		pcodethresh => 0.05
	);

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

	my $trigtime = $self->determine_trigger_time;

	my $stats=Util::HEAdas->new("ftstat") 
			      ->verbose(0);

	foreach my $file (@events) {

		# set TRIGTIME in EVENT extension

		if ($trigtime) {
			Util::FITSfile->new($file, 'EVENTS')
					->keyword(TRIGTIME => $trigtime,
						'[s] MET Trigger Time')
					;
		}
		
		my $evaux = 'NONE';
		$evaux = $filename->get('eventaux', 'b', '', 0)
		  if $file=~/bevshsp/;

		$tool->params({
					infile => $file,
					auxfile => $evaux,
				})
			->run;

		if($evaux ne 'NONE'){
		  my $auxfits = Util::FITSfile->new($evaux, 1);
		  $auxfits->specs('[col TIME]');
		  
		  $stats->params({infile  => $auxfits->fullname(),
				  centroid => 'no'})
		        ->run();

		  my $parfile = $stats->parfile();
		  my $start = Util::Date->new( $parfile->read('min') );
		  my $stop  = Util::Date->new( $parfile->read('max') );

		  $auxfits->ext(0);
		  $auxfits->specs('');

		  $auxfits->begin_many_keywords();
		  $auxfits->keyword('TSTART', $start->seconds());
		  $auxfits->keyword('TSTOP', $stop->seconds());
		  $auxfits->keyword('DATE-OBS', $start->date() .'T'.$start->time() );
		  $auxfits->keyword('DATE-END', $stop->date() .'T'.$stop->time() );
		  $auxfits->end_many_keywords();
		}
	}

}


1;