Subs::Filter (version 2.0)


package Subs::Filter;
##############################################################################
#
# DESCRIPTION: Filter the event data. First calculate attitude and
# DESCRIPTION: orbit related filtering quantities using the prefilter
# DESCRIPTION: tool. Then create the filter file by collating various
# DESCRIPTION: sources of HK with the prefilter output.
# DESCRIPTION: Next apply filtering criteria - spacial, temporal, or
# DESCRIPTION: event selection - to the event data.
# DESCRIPTION:
#
# HISTORY:
# HISTORY: $Log: Filter.pm,v $
# HISTORY: Revision 1.48  2006/09/10 20:05:16  apsop
# HISTORY: Remove explicit mode parameter setting from xrtscreen.
# HISTORY:
# HISTORY: Revision 1.47  2006/08/09 14:53:32  apsop
# HISTORY: Sort the xrt header file before using in making the filter file.
# HISTORY:
# HISTORY: Revision 1.46  2006/08/08 14:52:47  apsop
# HISTORY: Fix bugs in sorting of the DAP files.
# HISTORY:
# HISTORY: Revision 1.45  2006/08/04 13:36:00  apsop
# HISTORY: Check for duplicate time rows in dap files and remove.
# HISTORY:
# HISTORY: Revision 1.44  2006/05/10 15:20:13  apsop
# HISTORY: Remove RAWXTL column from filtered event lists.
# HISTORY:
# HISTORY: Revision 1.43  2006/04/27 16:21:59  apsop
# HISTORY: Comment out leapsec and tprec parameters from newmakefilter call. Allow for either INDEF or NULL values returned from FITSfile class.
# HISTORY:
# HISTORY: Revision 1.42  2006/01/31 16:44:48  apsop
# HISTORY: Switch to using snapshot gti to filter mkf file.
# HISTORY:
# HISTORY: Revision 1.41  2005/12/02 15:42:09  apsop
# HISTORY: Going back to using newmakefilter.
# HISTORY:
# HISTORY: Revision 1.40  2005/12/01 20:54:20  apsop
# HISTORY: Replace newmakefilter task with makefilter task.
# HISTORY:
# HISTORY: Revision 1.39  2005/11/20 20:31:18  apsop
# HISTORY: Check for presence of ACS_DATA extension in attitude file.
# HISTORY:
# HISTORY: Revision 1.38  2005/11/17 16:57:39  apsop
# HISTORY: Fix another typo, and clean up temporary files.
# HISTORY:
# HISTORY: Revision 1.37  2005/11/17 14:14:56  apsop
# HISTORY: Fix typo error.
# HISTORY:
# HISTORY: Revision 1.36  2005/11/08 17:38:57  apsop
# HISTORY: Fix bugs in gett s/c hk and bat dap hk files. Obtain min/max/median values for proper dap hk columns.  Make xrt event files of calibration only sources.
# HISTORY:
# HISTORY: Revision 1.35  2005/08/30 14:00:10  apsop
# HISTORY: Switch to using newmakefilter. Add proper NULL values for unsigned integer columns. Append used GTI to mkf file.
# HISTORY:
# HISTORY: Revision 1.34  2005/07/29 14:06:12  apsop
# HISTORY: remove bias map entries from xrt hd before using to make a filter file.
# HISTORY:
# HISTORY: Revision 1.33  2005/05/31 20:04:40  apsop
# HISTORY: Perform GTI filtering on .mkf file.  Derive min/max time using stats
# HISTORY: instead of dumping entire file.
# HISTORY:
# HISTORY: Revision 1.32  2005/05/04 15:19:14  apsop
# HISTORY: Give up on old algorithm and test for columns in cleaned xrt event list before deleting.
# HISTORY:
# HISTORY: Revision 1.31  2005/05/03 15:43:04  apsop
# HISTORY: Change globs for xrt modes that have b? instead of w?.
# HISTORY:
# HISTORY: Revision 1.30  2005/04/29 15:41:26  apsop
# HISTORY: Large change in setting of xrtscreen parameters. uvotscreen disabled. Remove unneeded columns from screened xrt event lists.
# HISTORY:
# HISTORY: Revision 1.29  2005/04/06 15:45:03  apsop
# HISTORY: Change to using CALDB for cal parameters.
# HISTORY:
# HISTORY: Revision 1.28  2005/03/25 19:55:03  apsop
# HISTORY: Fix bug in merging of HK data into mkf file.
# HISTORY:
# HISTORY: Revision 1.27  2005/03/15 19:57:56  apsop
# HISTORY: Fix bug in call for getting eng hk file names.  Write DATE-* keywords to mkf primary header.
# HISTORY:
# HISTORY: Revision 1.27  2005/03/15 19:03:26  apsop
# HISTORY: Check for presence of xrt no position message when making tdrss cat file.
# HISTORY:
# HISTORY: Revision 1.26  2005/02/08 18:18:38  apsop
# HISTORY: Move calculation of attorb file to Filter class. Split making of attitude flag columns in to two parts.
# HISTORY:
# HISTORY: Revision 1.25  2005/01/24 18:58:11  apsop
# HISTORY:
# HISTORY: Add support for SAFEHOLD column to makefilter file.
# HISTORY:
# HISTORY: Revision 1.24  2005/01/12 17:21:02  apsop
# HISTORY: Incorporate removal of _PNT suffix from output columns of prefilter
# HISTORY:
# HISTORY: Revision 1.23  2004/11/19 21:46:48  apsop
# HISTORY: New version of xrt2fits.
# HISTORY:
# HISTORY: Revision 1.22  2004/10/12 16:19:49  apsop
# HISTORY: Turn on data filtering to produce cleaned event lists
# HISTORY:
# HISTORY: Revision 1.21  2004/08/30 13:21:51  apsop
# HISTORY: Bobs new version of module, but comment out actual screening for now.
# HISTORY:
# HISTORY: Revision 1.20  2004/08/27 18:38:13  apsop
# HISTORY: Modified to use Ed's code to determine the filtering expressions, but
# HISTORY: the instrument tools for the actual filtering.
# HISTORY:
# HISTORY: Revision 1.19  2004/05/06 20:02:34  dah
# HISTORY: Add version number back into the header comments.
# HISTORY:
# HISTORY: Revision 1.18  2004/04/16 20:21:18  dah
# HISTORY: Begin using embedded history records
# HISTORY:
#
# VERSION: 2.0
#
#
##############################################################################


use Subs::Sub;

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

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

    $self->{DESCRIPTION}="Cleaning and filtering the event files";

    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();

    ##################################################
    # plot interesting quantities in the filter file
    ##################################################
  #  $self->make_plots();

    ##################################################
    # screen out bad event data
    ##################################################
    $self->filter_data();


} # 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');
      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");
           }
       }
    }

} # end of make_filter_file method


##############################################################################
# Make a number of plots of interesting quantities in the filter file
##############################################################################
sub make_plots {
    my $self = shift;

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

    ########################
    # get the filter file
    ########################
    my $filter = $filename->get("filter", "s");
    unless( -f $filter) { return; }

    my $fits_filter = Util::FITSfile->new($filter, 1);

    $log->entry("Plotting interesting things in the filter file");

    ######################################
    # set up for plotting
    ######################################
    my %yparm=(attd=>"RA DEC ROLL ANG_DIST",
               weel=>"FWHEEL"              ,
               bat1=>"BDIHKGRBFOLLOW BDIHKSLEWING BDIHKSAA BDIHKVALID",
               bat2=>"BDIHKNGOODDETS BCAHKDOINGBLK BBRHKTRIGNUM",
               saaa=>"ORBIT_SAA ACS_SAA BDIHKSAA"     );

    my $command_tmp=$self->temp_file("plot_commands");

    my $fplot=Util::Ftool->new("fplot")
                         ->params({infile=>$filter,
                                   xparm   => "TIME",
                                   rows    => "-",
                                   device  => "/cps",
                                   pltcmd => "\@$command_tmp",
                                   offset  => "yes"});

    my $plot_output="pgplot.ps";


    my $mode;
    foreach $mode (keys %yparm) {
        #################################################
        # make sure the filter file has all the columns
        # that we need for this plot
        #################################################
        my $good=1;
        foreach my $col (split /\s+/, $yparm{$mode}) {
            unless( $fits_filter->find_column($col)) {
                $log->entry("$filter is missing $col");
                $good=0;
                last;
            }
        }
        
        unless($good) {
            $log->entry("Skipping $mode plot");
            next;
        }

        #######################################
        # copy the plot command file and 
        # insert the file name for "MKF"
        #######################################
        my $template = $self->proctop()."/lists/filter_plot.$mode";
        open IN, "<$template";
        open OUT, ">$command_tmp";
        while (<IN>) {

            s/MKF/$filter/;
            print OUT;
        }
        close(OUT);
        close(IN);

        ##################################
        # make the plot
        ##################################
        unlink $plot_output;
        $fplot->params({yparm=>$yparm{$mode}})
              ->run();

        ############################################
        # save the plot with the correct file name
        ############################################
        if( -s $plot_output ) {
            rename $plot_output,
                   $filename->get("filterplot","s",$mode,0);
        } else {
            $log->error(2,"$plot_output not created");
        }
    } # end of loop over plot types


} # end of make_plots method


###########################################################################
# apply standard filters to the event data
###########################################################################
sub filter_data {

    my $self = shift;

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

    $self->xrtscreen();
    $self->xrt_cal_events();

#######	$self->uvotscreen;

} # end of filter_data'

sub filter_dataEAP {

    my $self = shift;

    my $log     =$self->log();
    my $filename=$self->filename();
    my $procpar =$self->procpar();
    my $jobpar  =$self->jobpar();
    
    #######################################
    # find the filter file
    #######################################
    my $filter = $filename->get("filter", "s");
    if( ! -f $filter ) {
        $log->entry("No filter file, so can't do screening");
        return;
    }


    ###########################
    # fire up the extractor
    ###########################
    my $extractor = Util::Extractor->new();
    $extractor->verbose(2);

    ########################################
    # now filter the data
    # this is going to need drastic changes
    ########################################
    my $inst;
    foreach $inst ("bat", "xrt", "uvot") {

        $log->entry("Filtering $inst files");

        ################################################
        # determine the modes which should be filtered
        ################################################
        my @modes=();
        if(   $inst eq "xrt" ) { @modes=("pcw?po","wtw?po", "puw?po", "lrw?po") }
        elsif($inst eq "uvot") { @modes=("*")    }
        elsif($inst eq "bat" ) { @modes=("evshps", "evshsl", "evshpo")    }

        ########################
        # loop over modes
        ########################
        my $mode;
        foreach $mode (@modes) {
            #################################################
            # determine what kid of files we should use as
            # input for filtering
            ################################################
            my $unfiltered_type="unfiltered";
            if($inst eq "xrt" && substr($mode,0,2) ne "pc" ) {
                ########################################################
                # for the XRT use the reconstructed event files unless
                # this is photon counting  mode
                ########################################################
                $unfiltered_type = "reconst";
            }

            #########################################
            # get the unfiltered files
            #########################################
            my @unfs = $filename->get($unfiltered_type, $inst, $mode, '*');
            
            #############################################
            # check if there are any files in this mode
            #############################################
            if(@unfs) {
                $log->entry("Mode $mode");
            } else {
                $log->entry("No unfiltered files for mode $mode\n");
                next;
            }

            ###############################################
            # get the mode name understood by the xselect
            # mission database
            ###############################################
            my $mdb_mode;
            if($inst eq "bat" ) { $mdb_mode = "EVENT" }
            elsif($inst eq "uvot" ) { $mdb_mode = "EVENT" }
            elsif($inst eq "xrt" ) {
                my $mo = substr($mode,0,2);
                if(   $mo eq "pc" ) { $mdb_mode="PHOTON"   }
                elsif($mo eq "wt" ) { $mdb_mode="WINDOWED" }
                elsif($mo eq "pu" ) { $mdb_mode="PILEDUP"  }
                elsif($mo eq "lr" ) { $mdb_mode="LOWRATE"  }
            
            }

            ####################################
            # set the mode and criteria to use
            ####################################
            $extractor->instrument($inst, $mdb_mode);

            ########################################
            # set the time filtering criteria
            ########################################
            if($inst eq "bat" ) {
                ##################################
                # BAT - no screening criteria yet
                # eventually just SAA probably
                ##################################
                $extractor->criteria("BDIHKNGOODDETS>1 && ".
                                     "BDIHKSAA==0 && ".
                                     "BHVMon_Level_MAX>99", $filter);

            } elsif($inst eq "xrt") {
                ######################################
                # XRT
                ######################################
                my $expression = $self->get_xrt_filter_expression();
                $extractor->criteria($expression, $filter);

            } elsif($inst eq "uvot") {
                $extractor->criteria("ELV>10 && ORBIT_SAA==0", $filter);
            }

            #######################################################
            # if there's no good time then there's no point doing
            # any filtering for this instrument
            ######################################################
            unless($extractor->has_good_time()) {
                $log->entry("No GTIs so skipping all $inst $mode files");
                next;
            }

            ###################################
            # region filter to use
            ###################################
            my $region="";
            if($inst eq "xrt" && $mode =~ /^pc/ ) {
                ##############################################
                # For XRT photon counting mode we screen out
                # the calibration sources
                ##############################################
                $region = $filename->fetch_cal("regstd", $inst, "");
            }

            $extractor->region($region);

            
            ################################################
            # get the row filter to use, if any
            # This is a CFITSIO extended filename syntax
            # row filtering statement which which will be
            # appended to the filename
            ################################################
            my $row_selection="";
            if($inst eq "xrt") {
                #######
                # XRT
                #######
                $row_selection = $self->get_xrt_row_selection($mdb_mode);
            } elsif($inst eq "uvot") {
                #############
                # UVOT
                #############
                $row_selection="QUALITY == 0";
            }

            $extractor->row_selection($row_selection);


            #############################
            # loop over unfiltered files
            #############################
            foreach my $unf (@unfs) {

                ##################################################
                # get the name of the corresponding filtered file
                ##################################################
                my $evt = $filename->corresponding($unfiltered_type, 'event', $unf);

                $log->entry("Extracting $evt from $unf");

                ###########################
                # run the extractor
                ###########################
                $extractor->infiles($unf)
                          ->outfile($evt, "event")
                          ->run();

                ##########################################
                # some special error checking
                ##########################################
                if(! $extractor->had_error() && ! -f $evt) {
                    $log->error(2, "Filtered file $evt not created");
                    next;
                }

                if($extractor->had_error() ) {
                    if( -f $evt) {
                        $log->entry("Deleting $evt");
                        unlink $evt;
                    }
                    next;
                }

                ####################################################
                # make sure the screened file has some events in it
                ####################################################
                my $nevents = Util::FITSfile->new($evt,"EVENTS")
                                            ->keyword("NAXIS2");
                if($nevents) {
                    $log->entry("$nevents filtered events");
                } else {
                    $log->entry("Deleteing $evt since it has $nevents events");
                    unlink $evt;
                    next;
                }

            } # end of loop over files
        } # end of loop over modes
    } # end of loop over instruments

} # end of filter_data method

############################################################################
#
############################################################################
sub read_xrt_range_file {
    my $self = shift;
    my $range_file=shift;
    my @extensions = @_;

    my $log     =$self->log();
    my $filename=$self->filename();
    my $procpar =$self->procpar();
    my $jobpar  =$self->jobpar();
    
        ##############################################
    # loop over the extensions of the range file
    ##############################################
    my @expression=();
    foreach my $ext (@extensions) {

        my $fits = Util::FITSfile->new($range_file, $ext);

        $fits->cols("PARNAME", "RANGE");
        my %range=$fits->table();

        ########################################
        # now loop over all the parameters
        ########################################
        foreach my $param (keys %range) {
            ##########################################
            # get the range string and
            # apply any aliases for the column names
            ##########################################
            my $range = $range{$param};


            #####################################
            # write the expression for this row
            #####################################
            if($range =~ /,/ ) {
                #############################
                # we have a min and or a max
                #############################
                my ($min, $max) = split /\s*,\s*/, $range;
                if($min !~ /(INDEF|NULL)/ ) { push @expression, ("$param >= $min") }
                if($max !~ /(INDEF|NULL)/ ) { push @expression, ("$param <= $max") }
            } else {
                ###################################################
                # no comma means force the param to a single value
                ###################################################
                push @expression, ("$param == $range");
            }
        }


    }

    return  join(" && ", @expression);

} # end of read_xrt_range_file method
    


############################################################################
# create a filter expression from an XRT-style parameter range file.
# This method remembers the expression it returned the first time it was
# run and will return the same value on subsequent calls without wasting
# any effort.
############################################################################
sub get_xrt_filter_expression {
    my $self = shift;

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

    if($self->{XRT_EXPRESSION} ) { return $self->{XRT_EXPRESSION} }

    ##################################
    # get the range file
    ##################################
    my $range_file = $filename->fetch_cal('rangefile', 'xrt');
    
    my $expression = $self->read_xrt_range_file($range_file, "ATTRANGE", "INSTRANGE");
    
    #####################################
    # alias "SAA" to "ORBIT_SAA"
    #####################################
    $expression =~ s/\sSAA\s/ ORBIT_SAA /;

    $self->{XRT_EXPRESSION} = $expression;

    return $self->{XRT_EXPRESSION};


} # end of make_expression_from_range_file method


############################################################################
#
############################################################################
sub get_xrt_row_selection {

    my $self = shift;
    my $datamode = shift;

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



    ################################################
    # find the event filter range calibration file
    ################################################
    my $range_file = $filename->fetch_cal('graderange', 'xrt');
    
    ###########################################
    # determine which extension we should read
    ###########################################
    my $ext="";
    if($datamode eq "LOWRATE" || $datamode eq "PILEDUP" ) { $ext = "PDRANGE" }
    elsif(                       $datamode eq "PHOTON"  ) { $ext = "PCRANGE" }
    elsif(                       $datamode eq "WINDOWED") { $ext = "WTRANGE" }
    else { $log->error(2, "Unknown DATAMODE $datamode")  }
    
    ##########################################
    # read the range file
    ##########################################
    my $expression = $self->read_xrt_range_file($range_file, $ext);

    return $expression;

} # end of get_xrt_row_filter method



sub uvotscreen {

    my $self = shift;

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

    my $option = $jobpar->read('opt_uvot_screen');
    if ($option =~ /^n/) {
        $log->entry('UVOT event screening disabled');
        return;
    }

    my $attorb = $filename->get('attorb', 's');

    my $uscreen = Util::HEAdas->new('uvotscreen')
                     ->params({
                        attorbfile => $attorb,
                        aoexpr => 'ANG_DIST.lt.100',
                        evexpr => 'QUALITY==0',
                        badpixfile => 'CALDB'
                     })
           		->is_script(1);

    # iterate over UVOT unfiltered event files
    foreach my $unfFile ($filename->get('unfiltered', 'uvot', '*', '*')) {

        my $evtFile = $filename->corresponding('unfiltered', 'event', $unfFile);

        $log->entry("filtering $unfFile");

        $uscreen->params({
              infile => $unfFile,
              outfile => $evtFile,
              })
           ->run();
    }

} # end of uvotscreen


sub xrtscreen {
    my $self = shift;

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

    #######################################
    # find the filter file
    #######################################
    my $filter = $filename->get('filter', 's');
	if (not -f $filter) {
		$log->entry("No filter file, so can't do screening");
		return
	}

    my $xrtscreen = Util::HEAdas->new('xrtscreen')
				->params({
                                          mkffile => $filter,
					  gtiexpr => 'DEFAULT',
					  exprgrade => 'DEFAULT',
					  expr => 'DEFAULT',
					  clobber => 'no',
					  chatter => 3,
					  history => 'yes',
					  createattgti => 'yes',
					  createinstrgti => 'yes',
					  gtiscreen => 'yes',
					  evtscreen => 'yes',
					  obsmodescreen => 'yes',
					  gtifile => 'xrtscreen_gti.tmp',
					  usrgtifile => 'NONE',
					  hkrangefile => 'CALDB',
					  timecol => 'TIME',
					  outfile => 'DEFAULT',
					  gtiext => 'GTI',
					  evtrangefile => 'CALDB',
					  cleanup => 'yes'
				 })
				->is_script(1);

    ################################################
    # determine the modes which should be filtered
    ################################################
    my @modes= qw(pcw?po wtw?po pub?po lrb?po pcw?s[lt] wtw?s[lt] pub?s[lt] lrb?s[lt]);

    my $inst = 'xrt';
    ########################
    # loop over modes
    ########################
    foreach my $mode (@modes) {

        my $mo = substr($mode,0,2);

        #################################################
        # determine what kid of files we should use as
        # input for filtering
        ################################################
        my $unfiltered_type="unfiltered";
        if ($mo ne "pc") {
            ########################################################
            # for the XRT use the reconstructed event files unless
            # this is photon counting  mode
            ########################################################
            $unfiltered_type = "reconst";
        }

	if($mode =~ /s.lt./){
	  $xrtscreen->params({createattgti => 'no'});
	}else{
	  $xrtscreen->params({createattgti => 'yes'});
	}

        #########################################
        # get the unfiltered files
        #########################################
        my @unfs = $filename->get($unfiltered_type, $inst, $mode, '*');

        #############################################
        # check if there are any files in this mode
        #############################################
        if (@unfs) {
            $log->entry("Mode $mode");
        } else {
            $log->entry("No unfiltered files for mode $mode\n");
            next;
        }

        #############################
        # loop over unfiltered files
        #############################
        foreach my $unf (@unfs) {

            ##################################################
            # get the name of the corresponding filtered file
            ##################################################
            my $evt = $filename->corresponding($unfiltered_type, 'event', $unf);

            $log->entry("Extracting $evt from $unf");

            ###########################
            # run xrtscreen
            ###########################
	    $xrtscreen->params({
				infile => $unf,
				outfile => $evt,
			       })
	               ->run;

	    unlink(qw(xselect.hty defaults.def xrtscreen_gti.tmp));

	    if (not -f $evt) {
	      $log->entry("Missing $evt");
	      next;
	    }

            ####################################################
            # make sure the screened file has some events in it
            ####################################################
	    my $fits = Util::FITSfile->new($evt,"EVENTS");
            my $nevents = $fits->keyword("NAXIS2");
            if($nevents) {
                $log->entry("$nevents filtered events");

		#############################
		# Delete unnecessary columns
		#############################
		my @cols = qw(CCDFrame PHAS Amp PixsAbove ROTIME OFFSET RAWXTL);
		my $specs='[col ';
		foreach my $cname (@cols){
		  $specs .= "-$cname; " 
		    if $fits->find_column($cname)
		}
		$specs .= ']';

	        $fits->specs($specs);
		$fits->copy('xrt_events.tmp');
		unlink $evt;
		rename 'xrt_events.tmp', $evt;
            } else {
                $log->entry("Deleting $evt since it has $nevents events");
                unlink $evt;
                next;
            }

	  } # end foreach unfiltered

	} # end foreach mode

} # end of xrtscreen

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;
}


sub xrt_cal_events {
  my $self = shift;

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

  #####################################################
  # Select out the events from the calibration sources
  #####################################################
  my $spec = '[(STATUS & b1111000000000000) > b0]';
  foreach my $win ('w2', 'w3'){
    my @phot = $filename->get('unfiltered', 'x', "pc${win}*", '*');
    my $xrtcal = $filename->get('xrtcal', 'x', "pc${win}", 0);
    my @more_files;
    my $i=0;
    foreach my $file (@phot) {
      my $fits_file = Util::FITSfile->new($file, 'EVENTS', $spec);
      if( -f $xrtcal ){
        my $tmp = 'xrt_cal_' . $i++ . '.tmp';
        $fits_file->copy($tmp);
        push @more_files, $tmp .'[EVENTS]';
      }else{
        $fits_file->copy($xrtcal);
      }
    }

    ###############################################################
    # If there is more than one file per window type, then need to
    # merge them together and update timing keywords.
    ###############################################################
    if(@more_files){
      my $merge_file = 'xrtcal_merge.tmp';
      Util::HEAdas->new('ftmerge')
                  ->params({infile => $xrtcal .','. join(',', @more_files),
                            outfile => $merge_file})
                  ->run();

      my $del = Util::HEAdas->new('ftdelhdu')
                            ->params({outfile => 'none',
                                      confirm => 'yes',
                                      clobber => 'yes'});

      $del->params({infile => $merge_file .'[GTI]'})
          ->run();
      $del->params({infile => $merge_file .'[BADPIX]'})
          ->run();

      my $stats=Util::HEAdas->new("ftstat") 
                            ->params({infile  => $merge_file .'[EVENTS][col TIME]',
	                              centroid => 'no'})
		                ->verbose(0)
				->run();

      my $parfile = $stats->parfile();
      my $tmin = $parfile->read('min');
      my $tmax = $parfile->read('max');

      my $merge_fits = Util::FITSfile->new($merge_file, 'EVENTS');
      $merge_fits->keyword('TSTART', $tmin);
      $merge_fits->keyword('TSTOP', $tmax);

      $merge_fits->ext(0);
      $merge_fits->keyword('TSTART', $tmin);
      $merge_fits->keyword('TSTOP', $tmax);

      my $first = Util::Date->new($tmin);
      my $last  = Util::Date->new($tmax);
      $merge_fits->keyword('DATE-OBS', $first->date().'T'.$first->time() );
      $merge_fits->keyword('DATE-END', $last->date().'T'.$last->time() );

      unlink $xrtcal;
      rename $merge_file, $xrtcal;
    }

    unlink (grep s/^(xrt_cal.*\.tmp).*$/$1/, @more_files);

    if( -f $xrtcal){
      my $rows = Util::FITSfile->new($xrtcal, 'EVENTS')->nrows();
      unlink $xrtcal unless $rows;
    }
  }
}

1;