package Subs::Filter;
##############################################################################
#
# DESCRIPTION: Apply filtering criteria - spacial, temporal, or
# DESCRIPTION: event selection - to the event data.
# DESCRIPTION:
#
# HISTORY:
# HISTORY: $Log: Filter.pm,v $
# HISTORY: Revision 1.52 2007/11/08 17:13:35 apsop
# HISTORY: Do not try and get calibration events from an empty xrt event list.
# HISTORY:
# HISTORY: Revision 1.51 2007/09/27 17:26:01 apsop
# HISTORY: The column name PHAS0 should be PHASO, with an O rather than a zero.
# HISTORY:
# HISTORY: Revision 1.50 2007/09/25 20:10:43 apsop
# HISTORY: Remove PHAS0 column from the xrt cleaned event lists.
# HISTORY:
# HISTORY: Revision 1.49 2007/09/11 18:15:27 apsop
# HISTORY: Move production of filter file to a separate class. Clean up obsolete subroutines.
# HISTORY:
# 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();
##################################################
# plot interesting quantities in the filter file
##################################################
# $self->make_plots();
##################################################
# screen out bad event data
##################################################
$self->filter_data();
} # end of body 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 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 PHASO 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);
next unless $fits_file->nrows();
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;