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;