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.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) {
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
################################################
$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);
##########################################
# 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
}
########################################################
# 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);
if ( -f $framehead ) {
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 $framehead FRAME D D % $xrt_cols{$column}\n";
}
} # end if we have a frame header file
###################################################
# BAT DAP HK
###################################################
my $dap = $filename->get("hk", "b", "dp", 0);
if( -f $dap) {
################################################
# here are all the columns we are interested in
################################################
my %cols=(BHVCtrl_Level => "BAT High Voltage Control Levels",
BHVMon_Level => "BAT High Voltage Monitor Levels",
DM_HK_Temp => "BAT DM Side temperatures",
DM_HK_HVLkCurr => "BAT High Voltage Leakage Current");
###########################################################
# 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;";
foreach my $col (keys %cols) {
foreach my $type ("MIN", "MAX") {
my $newcol="${col}_${type}";
$spec .= " $newcol(1I)=$type($col);";
print CONFIG "$newcol $minmax DAP_HK D D % $cols{$col}\n";
}
}
$spec .= "]";
######################################################
# create a temporary file with the min/max values
#####################################################
my $fits = Util::FITSfile->new($dap, "DAP_HK");
$fits->specs($spec);
$log->entry("Extracting min/max values from $dap");
$fits->copy($minmax);
} # end if there is a BAT DAP HK file
###################################################
# ACS bit flags
###################################################
my $attitude = $filename->get("attitude", "s");
if(-f $attitude) {
##################################################
# 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");
my $acs_fits = Util::FITSfile->new($attitude, "ACS_DATA");
$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
################################################################
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("makefilter");
$log->entry("Creating filter file $filter");
$makefilter->params({configure => $config,
infileroot => "",
outfile => $filter,
leapsec => $filename->fetch_cal("leapsec"),
tprec => 0.001,
mission => "SWIFT",
clobber => "yes" });
$makefilter->run();
my $makefits = Util::FITSfile->new($filter);
my @time = $makefits->cols('TIME')->table();
my $first = Util::Date->new($time[0]);
my $last = Util::Date->new($time[-1]);
$makefits->ext(0);
# $makefits->keyword('TSTART', $time[0]);
# $makefits->keyword('TSTOP', $time[-1]);
# $makefits->keyword('TELAPSE', $time[-1]-$time[0]);
$makefits->keyword('DATE-OBS', $first->date().'T'.$first->time() );
$makefits->keyword('DATE-END', $last->date().'T'.$last->time() );
} # 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->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 ne "INDEF" ) { push @expression, ("$param >= $min") }
if($max ne "INDEF" ) { 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 $badpix = $filename->fetch_cal('badpix', 'u');
my $attorb = $filename->get('attorb', 's');
my $uscreen = Util::HEAdas->new('uvotscreen')
->params({
attorbfile => $attorb,
aoexpr => 'ANG_DIST.lt.100',
evexpr => 'QUALITY==0',
badpixfile => $badpix,
})
->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
}
###################################
# expression to use
###################################
my $expression = $self->get_xrt_filter_expression();
my $xrtscreen = Util::HEAdas->new('xrtscreen')
->params({
createinstrgti => 'yes',
hkrangefile => $filename->fetch_cal('rangefile', 'xrt'),
evtrangefile => $filename->fetch_cal('graderange', 'xrt'),
mkffile => $filter,
expr => $expression,
})
->is_script(1);
################################################
# determine the modes which should be filtered
################################################
my @modes= qw(pcw?po wtw?po puw?po lrw?po);
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";
}
{
my $mdb_mode = '';
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" }
my $gradexpr = $mdb_mode
? $self->get_xrt_row_selection($mdb_mode)
: 'none';
# createattgti is only yes in PD unit test
my $createattgti = $mo eq 'pd' ? 'yes' : 'no';
# need to save a GTI file for each mode since it is required
# for later xrtimage run
my $gtifile = $filename->get('gti', 'xrt', $mo, 0);
# evtscreen/gtiscreen is yes for all event modes, no for IM mode
my $screen = $mo eq 'im' ? 'no' : 'yes';
$xrtscreen->params({
exprgrade => $gradexpr,
createattgti => $createattgti,
gtifile => $gtifile,
gtiscreen => $screen,
evtscreen => $screen,
});
}
#########################################
# 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;
}
###################################
# region filter to use
###################################
my $region="";
if ($mo eq 'pc') {
##############################################
# For XRT photon counting mode we screen out
# the calibration sources
##############################################
$region = $filename->fetch_cal("regstd", $inst, "");
}
#############################
# 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));
if (not -f $evt) {
$log->entry("Missing $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("Deleting $evt since it has $nevents events");
unlink $evt;
next;
}
} # end foreach unfiltered
} # end foreach mode
} # end of xrtscreen
1;