Subs::StartEndTimes (version 0.0)


package Subs::StartEndTimes;
##############################################################################
#
# DESCRIPTION: This subroutine does a number of things invoving the time  
# DESCRIPTION: values in every FITS HDU for a list of file types.
# DESCRIPTION: 
# DESCRIPTION: If an extension had TSTART and TSTOP keywords, it checks to
# DESCRIPTION: be sure that neither is zero and that TSTART <= TSTOP.
# DESCRIPTION: It also sets the values for the DATE-OBS, TIME-OBS, DATE-END,
# DESCRIPTION: TIME-END keywords. Finally it keeps a running tally of the 
# DESCRIPTION: earliest TSTART and the latest TSTOP. These values are used
# DESCRIPTION: to mark start and finish of the entire observation.
# DESCRIPTION: 
# DESCRIPTION: Every TIME column is checked for whether it is in order.
# DESCRIPTION: The BAT bcdh files are not checked since the BAT_TELEM column
# DESCRIPTION: is not supposed to be in order.
# DESCRIPTION: 
# DESCRIPTION: Finally, for every GTI extension, the subroutine checks that 
# DESCRIPTION: no two GTIs overlap, that they are in order and that they
# DESCRIPTION: all have STOP > START.
#
# HISTORY: 
# HISTORY: $Log: StartEndTimes.pm,v $
# HISTORY: Revision 1.31  2007/01/31 20:20:31  apsop
# HISTORY: Fix bug causing start/end times to be used from files tagged as invalid
# HISTORY:
# HISTORY: Revision 1.30  2006/08/08 14:51:09  apsop
# HISTORY: Do not use shared repository files when calculating global start and end times.
# HISTORY:
# HISTORY: Revision 1.29  2006/06/28 19:07:49  apsop
# HISTORY: Fix bug in setting tstart and TSTOP in primary header of GTI files.
# HISTORY:
# HISTORY: Revision 1.28  2006/02/07 16:38:48  apsop
# HISTORY: Do not do anything with gzipped files.
# HISTORY:
# HISTORY: Revision 1.27  2005/11/15 22:26:08  apsop
# HISTORY: Do proper checking for overlapped GTIs. Update ONTIME for bat event lists EVENTS extension if needed.
# HISTORY:
# HISTORY: Revision 1.26  2005/11/08 19:11:21  apsop
# HISTORY: calculate start, stop, and on time values for GTI extensions.
# HISTORY:
# HISTORY: Revision 1.25  2005/04/22 15:38:36  apsop
# HISTORY: Test whether column has any rows before extracting and filtering time information.
# HISTORY:
# HISTORY: Revision 1.24  2005/04/19 15:27:14  apsop
# HISTORY: Fix bug which was causing some extensions to be skipped.
# HISTORY:
# HISTORY: Revision 1.23  2005/03/25 19:50:18  apsop
# HISTORY: Use info from other extensions as fallback for setting times in primary hdu.
# HISTORY:
# HISTORY: Revision 1.22  2005/03/15 19:01:17  apsop
# HISTORY: Process all fits files, instead of just selected types.  Remove code for setting UTCFINIT keyword.
# HISTORY:
# HISTORY: Revision 1.21  2005/03/07 20:58:51  apsop
# HISTORY: Check BAT files for times < 1000 and remove them.
# HISTORY:
# HISTORY: Revision 1.20  2004/10/12 16:28:09  apsop
# HISTORY: Turn off sort checking of TIME columns in order to decrease run time.
# HISTORY:
# HISTORY: Revision 1.19  2004/08/30 13:19:30  apsop
# HISTORY: Add in new trend types.  Need better way of doing this.
# HISTORY:
# HISTORY: Revision 1.18  2004/08/22 18:42:25  apsop
# HISTORY: Initial changes for new file classes and repository
# HISTORY:
# HISTORY: Revision 1.17  2004/07/12 13:45:40  apsop
# HISTORY: Add timedata type to check list
# HISTORY:
# HISTORY: Revision 1.16  2004/06/08 00:09:36  apsop
# HISTORY: Fix for handling case with no attitude info.
# HISTORY:
# HISTORY: Revision 1.15  2004/06/02 18:53:20  apsop
# HISTORY: Add uvot compression trend file to list of files to proccess
# HISTORY:
# HISTORY: Revision 1.14  2004/05/28 19:45:20  apsop
# HISTORY: Write CVSINIT keyword into all the fits files
# HISTORY:
# HISTORY: Revision 1.13  2004/05/06 20:02:34  dah
# HISTORY: Add version number back into the header comments.
# HISTORY:
# HISTORY: Revision 1.12  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}="Determining start and end times of the observation";

    return $self;
}

##################
# METHODS:
##################

sub body {
    my $self=shift;

    my $log     =$self->log();
    my $filename=$self->filename();
    my $fileinfo = $filename->{'INFO'};
    my $procpar =$self->procpar();
    my $jobpar  =$self->jobpar();

    my @types = ( keys %{$fileinfo->{'all'}} );

    my ($tstart, $tstop);

    ##########################################
    # get the maximum reasonable TSTART-TSTOP
    ##########################################
    my $max_duration = $procpar->read("max_duration");

    ########################
    # loop over file types 
    ########################
    my $type;
    foreach $type (@types) {

        $log->entry("Examining TSTART and TSTOP in all $type files");

        ###################################
        # loop over the files of this type 
        ###################################
        my $file;
        foreach $file ($filename->any($type) ) {
	    next if $file =~ /\.gz/;
	    my ($inst, $mode, $index) = $filename->parse($file, $type);
	    next if $fileinfo->{$inst}->{$type}->{'notfits'};

            my $fits = Util::FITSfile->new($file);

            ################################
            # loop over all HDUs
            ################################
            my $nhdus = $fits->nhdus();
            my $ext;
	    my ($fstart, $fstop) = (1E10, 0);
	    ###################################################
	    # Count backwards so we can use info in other hdus
	    # to set values in primary header
	    ###################################################
            for($ext=$nhdus-1; $ext>=0; $ext--) {

	        $log->entry("Examining $file\[$ext\]");
	        $fits->ext($ext);

                ###########################################
                # check if there are TSTART/TSTOP keywords
                ###########################################
                my $start = $fits->keyword("TSTART");
                my $stop  = $fits->keyword("TSTOP");

		if( $ext > 0 && $inst eq 'b' && $fits->find_column('TIME') && $fits->nrows() ){
		  $fits->cols('TIME');
		  $fits->rows('1');
		  my $time1 = $fits->table();
		  if( $time1 < 1000 ){
		    ######################################
		    # remove 'events' with times lt 1000
		    ######################################
		    $fits->specs('[TIME > 1000]');
		    $fits->copy('f1000.tmp');
		    rename 'f1000.tmp', $fits->name();
		    $fits->specs('');

		    $start = $fits->table();
		    $fits->keyword("TSTART", $start);
		    
		    $log->error([1, ZERO_TIME_REMOVED], 
			  'Events with times less than 1000 were removed from ' . $fits->fullname());
		  }
		  $fits->cols('-');
		  $fits->rows('-');
		}

                ###########################################
                # check if this is a GTI extension
                ###########################################
                if($ext > 0) {
                  #####################################################
                  # primary extension doesn't have an EXTNAME keyword
                  #####################################################
		  my $extname = $fits->keyword('EXTNAME');
		  unless(defined $extname){
                    ################################################
                    # sanity check - this can be taken out when the 
                    # FITS files are more stable
                    ################################################
		    $log->error(1, "$file has no EXTNAME in extension 1");

		  } else {
		    if($extname =~  /GTI/) {

		      ###########################################
		      # are the START columns in order?
		      ###########################################
		      unless($fits->cols('START')->isOrdered() ) {
			$log->entry("$file\[$ext\] START column is out of order. Will sort.");
			$fits->sort('START');
		      }

		      ###############################################
		      # are the individual GTIs of positive duration 
		      ###############################################
		      my ($overlap, $prev_stop) = (0, 0);
		      my %intervals = $fits->cols("START", "STOP")->table();
		      foreach( sort {$a <=> $b} (keys %intervals) ) {
			my ($start, $stop) = ($_, $intervals{$_});
                        if($stop <= $start ) {
			  $log->error(2, "Invalid GTI".
				      "START=$start STOP=$stop ".
				      "in $file\[$ext\]" );
                        }
		        if($start <= $prev_stop){
			  $log->entry("GTI START=$start STOP=$stop in $file\[$ext\] ".
				      "overlaps with previous row, STOP=$prev_stop" );
			  
			  $overlap = 1;
			}
			$prev_stop = $stop;
		      }

		      ###########################################
		      # are there overlaps
		      ############################################
		      if($overlap) {
			$log->entry("$file\[$ext\] Has overlapping GTIs. Will merge.");

			my $temp = 'merged_gti.tmp';
			Util::Ftool->new('mgtime')
			           ->params({ingtis => "$file\[$ext\] $file\[$ext\]",
					     outgti => $temp,
					     merge  => 'OR'})
				   ->run();

			Util::FITSfile->new($temp)
			              ->import_header("$file\[$ext\]", 'except', '')
			              ->append_to($file);
			
			
			Util::Ftool->new('fdelhdu')
			           ->params({infile => "$file\[$ext\]",
					     confirm => 'no',
					     proceed => 'yes'})
				   ->run();
			unlink $temp;
			%intervals = $fits->cols("START", "STOP")->table();

			$ext++;
		      }

		      ##################################
		      # Determine TSTART, TSTOP, ONTIME
		      ##################################
		      my $ontime = 0;
		      $start = 1E10 unless $start;
		      $stop = 0 unless $stop;
		      foreach (keys %intervals ) {
			$start = $_ if $_ < $start;
			$stop = $intervals{$_} if $intervals{$_} > $stop;
                        $ontime += $intervals{$_} - $_;
		      }

		      $fits->keyword('TSTART', $start);
		      $fits->keyword('TSTOP', $stop);
		      $fits->keyword('ONTIME', $ontime);
		      if($overlap && $inst eq 'b' && $type eq 'unfiltered'){
			my $curr_ext = $fits->keyword('EXTNAME');
			$fits->ext(1);
			if( $fits->keyword('EXTNAME') eq 'EVENTS' ){
			  $fits->keyword('ONTIME', $ontime);
			}
			$fits->ext($curr_ext);
		      }
		    }
		  }
                }

                if(defined $start && defined $stop ) {
		  ########################################
		  # make some sanity checks
		  ########################################
		  if($start == 0.0 || $stop == 0.0 ||
		     $stop < $start || $stop - $start > $max_duration ) {
		    $log->error(1, "Skipping $file\[$ext\] since it has ".
				"invalid TSTART=$start TSTOP=$stop");
		    next;
		  }

		  $log->entry("TSTART=$start TSTOP=$stop");

		  #######################################################
		  # keep a running tally of the files absolute start/end_times
		  #######################################################
		  if($start<$fstart) { $fstart = $start }
		  if($stop >$fstop ) { $fstop  = $stop  }
		  
		  $log->entry("File fstart=$fstart fstop=$fstop");
		  
		  #######################################
		  # write the DATE/TIME OBS/END keywords
		  #######################################
		  my $first = Util::Date->new($start);
		  my $last  = Util::Date->new($stop);
		  
		  $fits->keyword('DATE-OBS', $first->date().'T'.$first->time() );
		  $fits->keyword('DATE-END', $last->date().'T'.$last->time() );
		}elsif( $ext==0 && $nhdus > 1){
		  ##################################################
		  # Use info in other hdus to update primary header
		  ##################################################
		  my $first = Util::Date->new($fstart);
		  my $last  = Util::Date->new($fstop);
		  
		  $fits->keyword('TSTART', $fstart);
		  $fits->keyword('TSTOP', $fstop);
		  $fits->keyword('DATE-OBS', $first->date().'T'.$first->time() );
		  $fits->keyword('DATE-END', $last->date().'T'.$last->time() );		  
		}
                $log->entry("Checking $file\[$ext\]");

            } # end of loop over HDUs

	    #############################################################
	    # keep a running tally of the absolute start/end_times
	    # Don't include share respository files in global start/stop
	    #############################################################
	    unless( $fileinfo->{$inst}->{$type}->{'repository'} ){
	      if(!defined($tstart) || $fstart<$tstart) { $tstart = $fstart }
	      if(!defined($tstop ) || $fstop >$tstop ) { $tstop  = $fstop  }
	    }
		  
	    $log->entry("Observation tstart=$tstart tstop=$tstop");

        } # end of loop over files
    } # end of loop over file types
    
    #############################################################
    # check if we got start and stop times from the FITS files
    #############################################################
    unless(defined $tstart && defined $tstop) {
        ########################################################
        # nothing from the FITS files, so try getting
        # the start and stop times directly from the telemetry
        ########################################################
        $log->entry("Extracting observation start and end times from ".
                    "the telemetry secondary headers.");
                    
        my $squirt = Util::Tool->new($procpar->read("squirt"), "squirt");
        $squirt->stdin(1); #so that we don't redirect from /dev/null
        foreach my $telem ($filename->any("telemetry")) {
        
            $squirt->command_line("-e '\"%time2\\n\"; filter: 0;' < $telem 2>&1");
            $squirt->run();

            my @times=split /\s+/, $squirt->stdout();


            @times = sort { $a <=> $b } @times;

            my $max = $times[@times-1];
            if(! defined $tstop || $max > $tstop) { $tstop = $max}
            @times = grep {$_ > $tstop - $max_duration } @times;
            my $min = $times[0];

            $log->entry("In $telem min=$min max=$max");
            
            if(!defined $tstart || $min < $tstart) {$tstart = $min}

        }
    }


    ############################################
    # check if we got absolute start/stop times
    ############################################
    unless(defined $tstart) {
        $log->error(1, "Can't determine observation start time, ".
                       "setting to zero");
        $tstart=0.0;
    }

    unless(defined $tstop) {
        $log->error(1, "Can't determine observation end time, setting to zero");
        $tstop=0.0;
    }


    ####################################
    # convert mission time to date/time
    ####################################
    my $start_date = Util::Date->new($tstart);
    my  $stop_date = Util::Date->new($tstop );

    
    ######################
    # log the results
    ######################
    $log->entry("overall TSTART=$tstart = ".
                $start_date->date()."T". $start_date->time() );

    $log->entry("overall TSTOP=$tstop = ".
                $stop_date->date()."T". $stop_date->time() );

    #####################################
    # ... and record them in the job.par
    #####################################
    $jobpar->set({tstart  => $tstart,
                  tstop   => $tstop,
                  obsdate => $start_date->date(),
                  obstime => $start_date->time(),
                  enddate =>  $stop_date->date(),
                  endtime =>  $stop_date->time() });


} # end of body method