Subs::MakeFilter (version 1.0)


package Subs::MakeFilter;
##############################################################################
#
# DESCRIPTION: Create the filter file by collating various
# DESCRIPTION: sources of HK with the prefilter output.
# DESCRIPTION:
#
# HISTORY:
# HISTORY: $Log: MakeFilter.pm,v $
# HISTORY: Revision 1.3  2008/05/16 14:26:40  apsop
# HISTORY: Make good attitude GTI.
# HISTORY:
# HISTORY: Revision 1.2  2007/11/08 17:16:03  apsop
# HISTORY: Add STLOCKFL and STLOCKST columns to the filter file.
# HISTORY:
# HISTORY: Revision 1.1  2007/09/11 18:24:35  apsop
# HISTORY: New class for creating the filter file.  Split off from Filter.pm because needed sooner.
# HISTORY:
#
# VERSION: 1.0
#
#
##############################################################################


use Subs::Sub;

@ISA = ("Subs::Sub");
use strict;

sub new {
    my $proto=shift;
    my $self=$proto->SUPER::new();

    $self->{DESCRIPTION}="Create filter file";

    return $self;
}

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

sub body {
    my $self=shift;

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

    #########################
    # create the filter file
    #########################
    $self->make_filter_file();

} # end of body method


#############################################################################
# assemble everything together into a single filter file
#############################################################################
sub make_filter_file {
    my $self=shift;

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

    #############################################################
    # get the name of the filter file
    #############################################################
    my $filter = $filename->get("filter", "s");
    unlink $filter;
    
    $log->entry("Creating filter file $filter");

    ####################################
    # create the configuration file
    ####################################
    my $config = $self->temp_file("makefilter_config");
    unlink $config;
    open CONFIG, ">$config";

    #########################################################
    # collect attitude/orbit derived information
    #########################################################
    my $attorb = $filename->get("attorb", "s");
    if( -f $attorb) {
      $self->fixNulls($attorb, 'PREFILTER');

      print CONFIG "ELV $attorb PREFILTER D D % ".
	           "angle between Z axis and earth limb\n";

      print CONFIG "BR_EARTH $attorb PREFILTER D D % ".
                   "Angular distance from the sunlit earth\n";
                     
      print CONFIG "FOV_FLAG $attorb PREFILTER D D % ".
                   "0=sky; 1=dark earth; 2=bright earth\n";

      print CONFIG "SUNSHINE $attorb PREFILTER D D % ".
                   "is the sun shining on the satelite?\n";

      print CONFIG "ANG_DIST $attorb PREFILTER D D % ".
                   "degrees from nominal pointing\n";

      print CONFIG "RA $attorb PREFILTER D D % ".
                   "Right Ascension of satelite pointing\n";

      print CONFIG "DEC $attorb PREFILTER D D % ".
                   "Declination of satelite pointing\n";

      print CONFIG "ROLL $attorb PREFILTER D D % ".
                   "Roll angle of satelite pointing\n";

      print CONFIG "SAA $attorb PREFILTER D D % ".
                   "SAA derived from the orbit data\n";
                     
      print CONFIG "SAA_TIME $attorb PREFILTER D D % ".
                   "SAA time derived from the orbit data\n";
                     
      print CONFIG "COR_SAX $attorb PREFILTER D D % ".
                   "cutoff rigidity calculated as for the SAX\n";
                     
      print CONFIG "MCILWAIN_L $attorb PREFILTER D D % ".
                   "McIlwain L parameter for Earth's magnetic field\n";
                     
      print CONFIG "SUN_ANGLE $attorb PREFILTER D D % ".
                   "Angle between sun and pointing direction\n";

      print CONFIG "MOON_ANGLE $attorb PREFILTER D D % ".
                   "Angle between moon and pointing direction\n";

      print CONFIG "RAM_ANGLE $attorb PREFILTER D D % ".
                   "Angle between satelite motion and pointing\n";
    }

    #############################################################
    # collect information from the engineering HK
    #############################################################
    my $proto_config = $self->proctop()."/lists/engineering_hk_filters";
    open PROTO, "< $proto_config";
    
    my $last_inst="";
    my $hk;
    my @extnames;
    while(<PROTO>) {
        chomp;
        
        #####################
        # skip comments and blank lines
        #####################
        if(/^\s*#/ || /^\s*$/ ) { next }

        #################################
        # parse the current line
        #################################
        my ($inst, $mnemonic, $column, $apids, $comment ) = split /\s*\|\s*/, $_;

        if( $inst ne $last_inst) {
            ###############################
            # new instrument
            ###############################
            $last_inst = $inst;
            
            ###############################################
            # get the HK file for this instrument and
            # extract a list of the extensions it contains
            ################################################
	    if($inst eq 'swift'){
	      $hk = $filename->get('scenhk', $inst, '', 0);
	    }else{
	      $hk = $filename->get('enhk', $inst, '', 0);
	    }
	    $log->entry("Extracting filter quantites from $hk");

            @extnames=();
            if(-f $hk ) {
                @extnames = Util::FITSfile->new($hk)
                                        ->list_hdus();
            }
            
        } # end if this is a new instrument
        
        #######################################################
        # a mnemonic can be in multiple APIDS,
        # we need to determine which APIDS are present and
        # merge them if neccessary
        ######################################################
        $log->entry("Extracting $mnemonic, a.k.a. $column from APIDS $apids");
        my @apids = split /\s+/, $apids;
        my $list = Util::FITSlist->new();

        foreach my $apid (@apids) {
            my $ext = sprintf("hk%03xx001", $apid);

            if(grep {$_ eq $ext} @extnames ) {
                ##############################
                # we have data for this APID
                ##############################
                $list->add($hk, $ext);
            }
            
        }
        
        #################################################
        # if we don't have any data, then we will leave
        # this mnemonic out of the filter file.
        #################################################
        if($list->count() == 0 ) {
            $log->entry("No data for $mnemonic, a.k.a. $column");
            next;
        }
        
        #############################################
        # if we get here, then merge the data
        # don't worry, if we only have one APID, then
        # we just use that one extension
        #############################################
        my $tmp = $self->temp_file($mnemonic);
        my $merged=$list->merge($tmp, ' ', 'TIME', $mnemonic);

        ##########################################
        # sort the merged file if there was more
        # than one HK file
        ##########################################
        my $ext = $list->get_extension();
        if($merged eq $tmp) {
            ##########################
            # mulitple sources
            ##########################
            my $fits = Util::FITSfile->new($merged);
            $fits->sort("unique");
            $ext = 1; # merged data written to the first extension
        }

	$self->fixNulls($merged, $ext);
        ########################################################
        # write a row to the makefilter config file
        ########################################################
        print CONFIG "$mnemonic $merged $ext D D $column $comment\n";
    
    } # end of loop over mnemonics

    #######################################################
    # add XRT frame HK
    #######################################################
    my $framehead = $filename->get("hk", "xrt", "hd", 0);
    my $tmphead = 'framehead.tmp';
    
    if ( -f $framehead ) {
      $self->fixNulls($framehead, 1);
      Util::HEAdas->new('ftsort')
	          ->params({infile => $framehead . '[1][XRTMode!=9]',
			    outfile => $tmphead,
			    columns => 'TIME'})
                  ->run();

      my %xrt_cols=(CCDTemp => "CCD Temperature",
                    Vod1 => "output drain voltage for left amp 1",
                    Vod2 => "output drain voltage for left amp 2",
                    Vrd1 => "reference voltage for amp 1",
                    Vrd2 => "reference voltage for amp 2",
                    Vsub => "substrate bias voltage",
                    Vbackjun => "back junction bias voltage",
                    Baselin1 => "Baseline voltage for Signal Chain A",
                    Baselin2 => "Baseline voltage for Signal Chain B",
		    ENDTIME => "",
		    PixGtULD => "");


      #####################################
      # loop over all the columns
      #####################################
      foreach my $column (keys %xrt_cols) {
	print CONFIG "$column $tmphead FRAME D D % $xrt_cols{$column}\n";
      }
    } # end if we have a frame header file
    
    ###################################################
    # BAT DAP HK
    ###################################################
    my $dap = $filename->get("bdaphk", "b", "", 0);
    if( -f $dap) {
        my $fits = Util::FITSfile->new($dap, "DAP_HK");

        ################################################
        # here are all the columns we are interested in 
        ################################################
        my %cols=(DM_HK_Temp     => "BAT DM Side temperatures",
                  DM_HK_HVLkCurr => "BAT High Voltage Leakage Current",
		  DM_HK_Vthr     => "BAT XA1 Discriminator Threshold Voltage",
		  DM_Cnts_LLD    => "BAT LLD Count",
		  DM_Cnts_Evt    => "BAT Event Count",
		  DM_Cnts_MLD    => "BAT Multi channel hit Count",
		  BHVCtrl_Level  => "BAT High Voltage Control Levels",
                  BHVMon_Level   => "BAT High Voltage Monitor Levels");

	###########################################################################
	# Remove duplicate TIME columns. ftsort seems to allow sorting on a max of
	# 5 columns
	###########################################################################
	my $old_rows = $fits->nrows();

	Util::HEAdas->new('ftsort')
	            ->params({infile => $dap,
			      outfile => 'sorted.tmp',
			      columns => 'TIME,' . join(',', (keys %cols)[0,1,2,3])})
		    ->run();

	unlink $dap;
	Util::HEAdas->new('ftsort')
	            ->params({infile => 'sorted.tmp',
			      outfile => $dap,
			      columns => 'TIME',
			      unique => 'yes'})
		    ->run();

	my $new_rows = $fits->nrows();
	if($old_rows != $new_rows){
	  $log->error(1, "Duplicate time rows in DAP file $dap: " . ($old_rows-$new_rows) .
		      " rows removed.");
	}
	unlink 'sorted.tmp';
                  
        ###########################################################
        # These are vector columns and we are just interested in 
        # the min and max values, so generate a CFITSIO extended
        # filename syntax statement which will calculate these.
        # At the same time we put entries into the makefilter config file
        # for these.
        # Some day we might add a mean or median
        ##########################################################
        my $minmax = $self->temp_file("dap_hk_min_max_values");
        my $spec="[col TIME;";
	my $n=0;
        foreach my $col (keys %cols) {
            foreach my $type ('MIN', 'MAX', 'MEDIAN') {
                my $newcol="${col}_${type}";

                my $colspec .= " $newcol(1J)=$type($col);";
		if( length($spec)+length($colspec) > 255 ){
		  $spec .= "]";
		  $fits->specs($spec);
		  $log->entry("Extracting min/max values from $dap to $minmax");
		  $fits->copy($minmax);

		  $spec="[col TIME;";
		  $n++;
		  $minmax = $self->temp_file("dap_hk_min_max_values$n");
		}

		$spec .= $colspec;
                print CONFIG "$newcol $minmax DAP_HK D D % $cols{$col}\n";
            }
        }

        $spec .= "]";
        $fits->specs($spec);
	$log->entry("Extracting min/max values from $dap to $minmax");
        $fits->copy($minmax);
        

    } # end if there is a BAT DAP HK file

    
    ###################################################
    # ACS bit flags
    ###################################################
    my $attitude = $filename->get("attitude", "s");
    if(-f $attitude) {
        my $acs_fits = Util::FITSfile->new($attitude);
	if( (grep /ACS_DATA/, $acs_fits->list_hdus()) ){
	  $acs_fits->ext('ACS_DATA');

	  ##################################################
	  # Split the bit flags into separate columns in a
	  # temporary file
	  ##################################################
	  my $flags1 = $self->temp_file("acs_flags1");
	  my $flags2 = $self->temp_file("acs_flags2");

	  $acs_fits->specs("[col TIME; ".
			   "TEN_ARCMIN(1B)=FLAGS .eq. b1xxxxxxx ; ".
			   "SETTLED(1B)=FLAGS .eq. bx1xxxxxx]");

	  $acs_fits->extract($flags1);
        
	  $acs_fits->specs("[col TIME; ".
			   "ACS_SAA(1B)=FLAGS .eq. bxx1xxxxx ; ".
			   "SAFEHOLD(1B)=FLAGS .eq. bxxx1xxxx]");

	  $acs_fits->extract($flags2);
        
	  ################################################################
	  # now list the extracted columns to the config file
	  ################################################################
	  $self->fixNulls($flags1, 'ACS_DATA');
	  $self->fixNulls($flags2, 'ACS_DATA');

	  print CONFIG "TEN_ARCMIN $flags1 ACS_DATA D D % ACS reports within 10 arcmin of target\n";
	  print CONFIG    "SETTLED $flags1 ACS_DATA D D % ACS reports settled on target\n";
	  print CONFIG    "ACS_SAA $flags2 ACS_DATA D D % ACS reports in SAA\n";
	  print CONFIG   "SAFEHOLD $flags2 ACS_DATA D D % ACS reports in SAFE\n";
	}
    } # end if there is an attitude file


    #################################################
    # that's it for the config file, so close it up
    #################################################
    close CONFIG;

    ################################################
    # check if there is anything in the config file
    ################################################
    unless( -s $config ) {
        $log->entry("Nothing to put in the filter file");
        unlink $config;
        return;
    }


    ##############################
    # create the filter file
    ##############################
    my $makefilter = Util::HEAdas->new("newmakefilter");
    $log->entry("Creating filter file $filter");
    $makefilter->params({configure  => $config,
                         infileroot => "",
                         outfile    => $filter,
##                         leapsec    => 'CALDB',
##                         tprec      => 0.001,
                         mission    => "SWIFT",
                         clobber    => "yes"   });

    $makefilter->run();
    unlink $tmphead if -f $tmphead;

    {
      my $makefits = Util::FITSfile->new($filter, 'FILTER');
      #
      # Calculate additional star tracker info
      # 
      $makefits->specs('[col STLOCKFL(1B)=(SAC_MODESTAT/32)%2==1; '.
		       'STLOCKST(1I)=STA0VALID+STA1VALID+STA2VALID+STA3VALID+STA4VALID+STA5VALID; '.
		       '-STA0VALID; -STA1VALID; -STA2VALID; -STA3VALID; -STA4VALID; -STA5VALID;]');
      $makefits->copy($tmphead);
      $makefits->specs('');
      rename $tmphead, $filter;

      $makefits->cols('STLOCKFL');
      $makefits->column_comment('Star Tracker locked flag');
      $makefits->cols('STLOCKST');
      $makefits->column_comment('Number of currently valid Star Tracker stars');

      my ($sum,$mean,$sigma,$min,$max) = $makefits->cols('TIME')->stats;
      $makefits->keyword('TSTART', $min);
      $makefits->keyword('TSTOP', $max);

      my $first = Util::Date->new($min);
      my $last  = Util::Date->new($max);

      $makefits->ext(0);
      $makefits->keyword('DATE-OBS', $first->date().'T'.$first->time() );
      $makefits->keyword('DATE-END', $last->date().'T'.$last->time() );
    }

    {
       my $gtis = $filename->get('gti', 's', 'ss');
       if ($gtis) {
           my $tmpfile = 'constrained.mkf';
           Util::HEAdas->new('ftcopy')
                  ->params({
                      infile => $filter . qq([1][gtifilter(\\"$gtis\\")]),
                      outfile => $tmpfile,
                  })->run;
           if (-f $tmpfile) {
               unlink($filter);
               rename($tmpfile, $filter)
                 or $log->error(1, "unable to rename $tmpfile to $filter [$!]");
	       Util::HEAdas->new('ftappend')
		   ->params({infile => "$gtis\[GTI\]",
			     outfile => $filter})
		   ->run();
           }
           else {
              $log->error(1, "error applying GTIs $gtis to $filter");
           }
       }

     }


    #########################
    # Make good attitude gti
    #########################
    my $gti = $filename->get('gti', 's', 'at');
    my $expr = 'SAFEHOLD==0 && SAC_ADERR<0.2 && STLOCKFL==1 && STAST_LOSSFCN<' . 
               $procpar->read('att_loss_thresh');

    my $maketime = Util::Ftool->new('maketime')
      ->params({infile => $filter.'[FILTER]',
		outfile => $gti,
		expr => $expr,
		prefr => 0,
		postfr => 0,
		compact => 'no'})
      ->run();

    my $gtifits = Util::FITSfile->new($gti);
    for(my $ext=0; $ext < 2; $ext++){
      $gtifits->ext($ext);
      $gtifits->keyword('MJDREFI', 51910, 'MJD reference day');
      $gtifits->keyword('MJDREFF', 7.428703700000000E-04, 'MJD reference (fraction of day)');
    }
				      
} # end of make_filter_file method


sub fixNulls {
  ############################################################
  # Set TNULL values for columns with unsigned integer types.
  # Need this as makefilter uses an inappropriate default 
  # TNULL = -max for these columns
  ############################################################
  my ($self, $file, $ext) = @_[0,1,2];

  my $fits = Util::FITSfile->new($file, $ext);
  my %keywords = $fits->keywords();
  my %max = ('B' => 256,
	     'I' => 32768,
	     'J' => 2147483648);

  $fits->begin_many_keywords();
  my $n=0;
  foreach my $tform (grep /^TFORM/, keys %keywords){
    next unless $keywords{$tform} =~ /([IJ])/;
    my $type = $1;
    my $num = $tform;
    $num =~ s/TFORM//;

    next unless $keywords{"TZERO$num"} = $max{$type};

    unless( $keywords{"TNULL$num"} ){
      $fits->keyword("TNULL$num", $max{$type}-1);
      $n++;
    }
  }
  $fits->end_many_keywords() if $n;
}

1;