package Util::Extractor;
##############################################################################
#
# DESCRIPTION: This is an interface to the extractor FTOOL which is the standard
# DESCRIPTION: program for screening event files and extracting images,
# DESCRIPTION: spectra and lightcurves.
# DESCRIPTION: This class does a little more than just run the extractor.
# DESCRIPTION: It can also handle GTI files and knows mission specific
# DESCRIPTION: information by reading the xselect.mdb file.
# DESCRIPTION: This class should be able to replace calls to xselect and
# DESCRIPTION: do so more efficiently and with better error checking.
#
# HISTORY
# HISTORY: $Log: Extractor.pm,v $
# HISTORY: Revision 1.2 2006/08/01 20:35:34 apsop
# HISTORY: Add in CVS history indicator.
# HISTORY:
# HISTORY: 1.0 -> 1.1 2000-06-20
# HISTORY: Now the instrument method returns $self.
# HISTORY: Fixed a bug in the gti method which mishandles a single gti file name
# HISTORY:
# HISTORY: 1.1 -> 1.2 2000-08-09
# HISTORY: Now set binf to handle ASCA SIS region selection properly.
# HISTORY: Added bin_size method to set lightcurve bin size.
# HISTORY: Made temporary file names unique.
# HISTORY:
# HISTORY: 1.2 -> 1.3 2002-04-22
# HISTORY: Slight modification of initialization to make it work with
# HISTORY: new Ftools class
# HISTORY:
# HISTORY: 1.3 -> 1.4 2003-05-28
# HISTORY: Add support for 'raw' coordinate type
# HISTORY:
# HISTORY: 1.4 -> 1.5 2003-07-10
# HISTORY: Now correctly handle the situation where the screening
# HISTORY: criteria produice no GTIs. This required a change to the
# HISTORY: gti method. Now you specify "NONE" to that method to
# HISTORY: indicate that all times are good. No arguments indicates
# HISTORY: that no times are good. Added the has_good_time method.
# HISTORY: The extrator will not run if there is no good time.
# HISTORY:
# HISTORY: 1.5 -> 1.6 2004-04-08
# HISTORY: now read the "events" parameter from the MDB
# HISTORY:
# HISTORY: 1.6 -> 1.7 2004-05-26
# HISTORY: now explicitly clode the MDB file after reading it
# HISTORY:
# HISTORY: 1.7 -> 1.8 2004-06-04
# HISTORY: now set the ccol parameter
#
# VERSION: 1.8
#
##############################################################################
use Util::Ftool;
@ISA=("Util::Ftool");
use strict;
use Util::FITSfile;
#################################
# used to keep file names unique
#################################
my $INVOCATION;
###############################
# mission database class data
###############################
my $MDB={}; # mission database
my $KEYS=""; # a list of unique keys in the mission database
my %ALIASES=(ae=>"ASTRO-E",
xrs=>"XRS",
xis0=>"XIS-0",
xis1=>"XIS-1",
xis2=>"XIS-2",
xis3=>"XIS-3",
hxd_wel=>"HXD-WELL",
hxd_trn=>"HXD-TRN",
mn=>"ASCA",
ad=>"ASCA",
s0=>"SIS0",
s1=>"SIS1",
g2=>"GIS2",
g3=>"GIS3",
sw=>"SWIFT",
st=>"SWIFT",
xrt=>"XRT",
uvot=>"UVOT",
bat=>"BAT" );
#there's probably a better way of doing this
################################################
# constructor
################################################
sub new {
my $self=shift;
$self=$self->SUPER::new("extractor");
##################################################################
# increment the invocation count and set the temp file extension
##################################################################
$INVOCATION++;
my $tmp_ext="_extractor_${INVOCATION}.tmp";
###################################
# temporary file names
###################################
$self->{REGION_TMP}="region$tmp_ext";
$self->region();
$self->{GTI_TMP}="gti$tmp_ext";
$self->{NO_GOOD_TIME}=0;
$self->gti("NONE");
$self->{ROW_SELECTION}="";
$self->{ROW_SELECTION_FILE}="";
$self->{INFILES_TMP}="infiles$tmp_ext";
unlink $self->{INFILES_TMP};
##################################################
# flag to set if any of the criteria are
# so invalid that we should not run the extractor
##################################################
######################################
# initially empty instrument database
######################################
$self->{IDB}={};
##############################################
# maximum length of a criterion string before
# it must be written to a temporary file
##############################################
$self->{CRITERIA_INLINE_LENGTH}=80;
$self->{GTIS_INLINE_LENGTH}=80;
############################
# default coordinate types
############################
$self->{COORD_TYPE} = "det";
$self->{WMAP_COORD_TYPE} = "det";
###########################################
# name of a sample input file which
# can be queried if we need to know about
# e.g. column order
###########################################
$self->{SAMPLE_INPUT}="";
return($self);
}
#####################
# ACCESSORS:
#####################
#####################################################################
# returns false if all the time has been filtered out.
# returns true if at least some good time made it through the filter.
#####################################################################
sub has_good_time {
my $self = shift;
return ! $self->{NO_GOOD_TIME};
}
########################################################################
########################################################################
# if the Util::Tool verbosity is greater than one and the log is defined,
# returns the log object. Otherwise returns false
########################################################################
sub very_verbose {
my $self=shift;
my $log=$self->log();
if($log && $self->verbose_level()>1) { return $log }
else {return undef }
}
########################################################################
########################################################################
# Set the mission parameters in the extractor class data.
# This method reads and parses the xselect.mdb mission database file.
# Note you have to also set the instrument before you can apply
# any of these values.
########################################################################
sub mission {
my $self=shift;
#####################################
# determine mission alias if any
#####################################
my $chosen_mission=$self->alias(shift);
###########################################################
# empty the mission database if it has already been filled
###########################################################
if($MDB) { $MDB={}; };
################################################
# find the xselect mission database file
################################################
my $mdb_file=$self->install_dir()."/bin/xselect.mdb";
unless(-f $mdb_file ) {
$self->setup_error(3, "No mission databse file $mdb_file" );
}
################################################
# open the file containing the mission database
# and loop over all its rows
################################################
open MDB, "<$mdb_file" or return;
my $tag;
my $key;
my $value;
my $mission;
my %unique_keys;
while(<MDB>) {
chomp;
#################################
# ignore blank and comment lines
#################################
unless( /^\s*!/ || /^\s*$/ ) {
###############################
# parse the line
###############################
($tag,$value) = /^(\S*)\s*(.*)$/;
my @inst=split /:/, $tag;
$mission=shift @inst;
$key= pop @inst;
my $inst=join ":", @inst;
if($mission eq $chosen_mission ) {
####################################################
# initialize the MDB entry for this key if it does
# not already exist
####################################################
unless($MDB->{$key}) { $MDB->{$key}={} }
###########################################
# save this keyword in the mission database
###########################################
$MDB->{$key}->{$inst}=$value;
###########################################
# update the list of unique keys
###########################################
$unique_keys{$key}="";
} # end if this line if for our chosen mission
} # end of this is not a comment or blank line
} # end of loop over lines in xselect.mdb file
close MDB;
###################################
# save the list of unique keywords
###################################
$KEYS = [ keys %unique_keys];
} # end of mission method
##########################################################################
# Set the instrument-specific parameters for this extractor object.
# It actually populates the instrument database by extracting values
# from the mission database.
##########################################################################
sub instrument { #(instrument, mode, submode ...)
my $self=shift;
#####################################################
# assemble all mission specifiers into a single array
# and apply aliases
#####################################################
my @instrument=split /:/, join(":", (@_));
foreach (@instrument) {
$_ = $self->alias($_);
}
#####################################################
# loop over all the keywords in the mission database
#####################################################
my $key;
foreach $key ( @$KEYS ) {
#################################################################
# loop over sucessively more general instrument specifications
# until we get to one which has an entry.
#################################################################
my $value;
my @inst=@instrument;
while(! defined ($value=$MDB->{$key}->{join ":", @inst}) ) {
if( !@inst) {
#############################################
# there are no matches for this keyword,
# so break out of the loop.
# Otherwise, we would loop forever.
##############################################
last;
}
pop @inst;
} # end of loop over instrument specifier generalizations
#####################################################
# save the keyword values in the instrument database
#####################################################
if(defined $value) {
$self->{IDB}->{$key} = $value;
}
} # end of loop over keywords
return $self;
} # end of instrument method
#############################################################################
# get a value out of the mission database. Actually this method gets a value
# from the instument databse which is a subset of the mission database.
# The mission and instrument must be set before this method may be used
#############################################################################
sub mdb {
my $self=shift;
my $key=shift;
return $self->{IDB}->{$key} || "";
} # end of mdb method
###########################################################################
# get an alias for a mission or instrument name.
# This translates from the terms we will commonly use to the terms
# used in the xselect.mdb file
###########################################################################
sub alias {
my $self=shift;
my $value=shift;
if(exists $ALIASES{$value}) { return $ALIASES{$value} }
else { return $value }
} # end of alias method
##########################################################
# Set the region filters to be used in the extractor run.
# Concatenates multiple files.
##########################################################
sub region { #(region1, region2 ...)
my $self=shift;
my @list=@_;
################################
# clear out any temporary files
################################
unlink $self->{REGION_TMP};
if(@list==0 || !$list[0] ) {
####################
# no region filter
####################
$self->params({regionfile=>"NONE"});
return $self;
}
my $regionfile;
if(@list==1) {
################################################
# single region filter - specify it explicitly
###############################################
$regionfile=$list[0];
} else {
################################################
# multiple regions. We have to concatenate them
################################################
$regionfile=$self->{REGION_TMP};
open OUT, "$regionfile";
foreach (@list) {
if( ! -e $_ ) {
##########################################
# file does not exist, so report an error
##########################################
$self->setup_error(2,"Region filter $_ does not exist");
} else {
#####################################
# concatenate file
#####################################
open IN, "<$_";
print OUT <IN>;
close IN;
}
} # end of loop over files
close OUT;
}
##################################
# set the parameter
# and talk to the log
##################################
$self->params({regionfile=>$regionfile});
my $log;
if($log=$self->very_verbose() ) {
$log->entry("Filtering with region file:");
$log->file($regionfile);
}
return $self;
}
########################################################################
# set the GTI time screening files.
# Multiple GTIs are combined with the specified logic ("and" or "or").
# specify "NONE" to indicate that no time filtering should be done.
# specify no arguments to indicate that there is no good time
# and that the filtering should produce no output.
#######################################################################
sub gti { #(gti1, gti2, logic ...)
my $self=shift;
#############################################
# reset the "all time filtered out" flag
#############################################
$self->{NO_GOOD_TIME}=0;
##############################################
# no arguments means no good time
# we need to set a flag to indicate that we will
# not be filtering
##############################################
if( @_ == 0 ) {
$self->{NO_GOOD_TIME}=1;
$self->params({timefile=>""});
return $self;
}
##############################################
# "NONE" means don't time filter
##############################################
if( $_[0] eq "NONE" ) {
$self->params({timefile=>"NONE"});
return $self;
}
###############################################
# If there are multiple GTIs given, the last
# argumement should be the logic indicator
# it may be indicated with several formats and
# defaults to "OR"
#############################################
my $logic=pop;
if( uc($logic) eq 'AND' || $logic eq '&&' || $logic eq '&') {
################
# AND
################
$logic="AND";
} elsif(uc($logic) eq 'OR' || $logic eq '||' || $logic eq '|') {
###############
# OR
###############
$logic="OR";
} else {
################################################
# this looks more like a GTI file name
# so set to the default and push the argument
# back onto the GTI list
################################################
push @_, ($logic);
$logic="OR";
}
###################################################
# the arguments left in @_ are the GTI file names
# clear out any temporary files
###################################################
my @gtis=@_;
my $timefile;
if($gtis[0] ne $self->{GTI_TMP}) {unlink $self->{GTI_TMP} };
if(@gtis==1) {
################################################
# single GTI file - specify it explicitly
################################################
$timefile=$gtis[0];
} else {
################################################
# multiple GTI files. We have to merge them
################################################
my $ingtis = join " ", @gtis;
my $tmp_file;
if(length($ingtis) > $self->{GTIS_INLINE_LENGTH} ) {
#########################################################
# we need to dump the GTIs to a file for input to mgtime
#########################################################
$tmp_file="extractor_gti_list.tmp";
open GTIS, ">$tmp_file";
print GTIS join "\n",@gtis;
close GTIS;
$ingtis="\@$tmp_file";
}
####################################
# merge
####################################
$timefile=$self->{GTI_TMP};
Util::Ftool->new("mgtime")
->params({ingtis => $ingtis,
outgti => $timefile,
merge => $logic,
instarts => "START",
instops => "STOP",
indates => "MJDREF",
intimes =>" ",
outstart => "START",
outstop => "STOP"})
->run();
#################################
# cleanup
#################################
if( defined $tmp_file) {unlink $tmp_file}
}
#######################################################
# set the parameter to point to the combined GTI file
# set the GTI extension name to STDGTI
# since this is what maketime outputs
# There's probably a more general and
# elegant way of handling this
#######################################################
$self->params({timefile => $timefile,
gtinam => "STDGTI" });
####################################
# dump the GTIs to the log
####################################
my $log;
if($log=$self->very_verbose() ) {
$log->entry("GTIs from screening criteria:");
$log->text(Util::FITSfile->new($timefile)
->dump_table() );
}
return $self;
}
#############################################################
# set the GTIs by way of a set of screening criteria
# and a filter file by running the maketime FTOOL
#############################################################
sub criteria {
my $self=shift;
my $criteria=shift;
################################################
# log the screening criteria
################################################
my $log;
if($log=$self->very_verbose() ) {
$log->entry("Using screening criteria:");
$log->text($criteria);
}
###############################
# set up the maketime FTOOL
###############################
my $maketime=Util::Ftool->new("maketime")
->params({compact=> "no",
copykw => "yes",
histkw => "no",
prefr => 0.0,
postfr => 1.0});
##################################################################
# We will place short criteria in the parfile directly
# but long or multiple line ones have to go into a temporary file
##################################################################
my $criteria_file="extractor_criteria.tmp";
if( $criteria !~ /\n/ &&
length($criteria) <=$self->{CRITERIA_INLINE_LENGTH} ) {
########################
# give criteria in-line
########################
$maketime->params({expr=>$criteria});
} else {
######################################
# write criteria to a temporary file
######################################
unlink $criteria_file;
open CRITERIA, ">$criteria_file";
print CRITERIA "$criteria";
close CRITERIA;
$maketime->params({expr=>"\@$criteria_file"});
}
################################################################
# loop over all filter files and make a GTI file for each
################################################################
my @gtis=();
my $filter;
while( $filter=shift ) {
##########################################
# run maketime
##########################################
my $gti="${filter}_gti.tmp";
unlink $gti;
$maketime->params({infile => $filter,
outfile => $gti })
->run();
###################################################
# did we produce a GTI file with something in it?
###################################################
if(! $maketime->had_error() && Util::FITSfile->new($gti,1)
->keyword("NAXIS2") ){
################################################
# this is a good GTI file, so remember its name
################################################
push @gtis, ($gti);
} else {
###############################
# no GTIS in this file
###############################
unlink $gti;
my $log;
if($log=$self->very_verbose() ) {
$log->entry("No GTIs from $filter");
}
} # end of check for whether we produced a good GTI file
} # end of loop creating a GTI for each filter
###############################################################
# check how many GTI files we just made and behave accordingly
###############################################################
if(@gtis == 0 ) {
################################################
# no GTI files
################################################
$log->entry("No GTIs produced");
$self->gti();
} elsif( @gtis == 1) {
#####################################
# there's only one GTI, so rename it
# to the official temp file
####################################
rename $gtis[0], $self->{GTI_TMP};
$gtis[0] = $self->{GTI_TMP};
$self->gti(@gtis);
} else {
###################################################
# multiple GTIs - need to clean up temporary files
###################################################
$self->gti(@gtis,"OR");
foreach (@gtis) {unlink;}
}
#######################
# cleanup
#######################
unlink $criteria_file;
return $self;
} # end of criteria method
#################################################################################
# This method sets a row selection criterion to apply to all the input
# files. It should be in the form of a CFITSIO extended filename
# row selection statement, except without enclosing brackets.
# The expression may contain newline characters.
# If the argument is blank or omitted, then any previously set criteria
# will be removed.
# This method must be called before setting the input files.
################################################################################
sub row_selection {
my $self = shift;
my $criteria=shift || "";
#################################
# report the criteria to the log
#################################
my $log;
if($criteria && ($log=$self->very_verbose()) ) {
$log->entry("Using row selection criteria:");
$log->text($criteria);
}
######################################################
# if there is no argument or the argument is blank,
# then we have no row filtering criteria
#####################################################
if(!$criteria) {
$self->{ROW_SELECTION} ="";
} elsif(length($criteria) > $self->{CRITERIA_INLINE_LENGTH} ||
$criteria =~ /\n/ ) {
#########################################
# long or multi-line criteria
# so dump them to a temporary file
#########################################
my $tmp = "extractor_row_selection_${INVOCATION}.tmp";
unlink $tmp;
open TMP, ">$tmp";
print TMP $criteria;
close TMP;
$self->{ROW_SELECTION} = "[ @".$tmp." ]";
$self->{ROW_SELECTION_FILE} = $tmp;
} else {
######################################################
# the row selection is simple enough we can inline it
######################################################
$self->{ROW_SELECTION} = "[ $criteria ]";
}
return $self;
} # end of row_selection method
####################################################
# input event file(s)
####################################################
sub infiles { #(file1, file2 ...)
my $self=shift;
my @list=@_;
if(@list == 0 ) {
#####################################
# no input files - this is an error
####################################
$self->setup_error(2,"Specified nothing for extractor input file");
} elsif(@list==1) {
##################################
# one input file list it directly
##################################
$self->params({filename=>$list[0].$self->{ROW_SELECTION}});
$self->{SAMPLE_INPUT} = $list[0];
} else {
################################################
# multiple input files - these must be listed
# in a temporary file
###############################################
open LIST, ">$self->{INFILES_TMP}";
foreach (@list) {
print LIST "${_}$self->{ROW_SELECTION}\n";
}
close LIST;
$self->params({filename=>"\@$self->{INFILES_TMP}" });
$self->{SAMPLE_INPUT} = $list[0];
}
return $self;
} # end of infiles method
#############################################################################
# set the lightcurve bin size and optionally the fraction of a bin
# which has to be within a GTI for it to be included. This defaults
# to "1".
#############################################################################
sub bin_size { #(bin_size, lcthresh )
my $self=shift;
my $bin_size=shift;
my $lcthresh=shift;
unless(defined($lcthresh)) {$lcthresh=1}
$self->params({binlc => $bin_size,
lcthresh => $lcthresh });
return $self;
} # end of bin_size method
##############################################################################
# set the coordinates "det" or "sky" to be used for filtering and
# image extraction. The default is "det".
##############################################################################
sub coord_type {
my $self=shift;
my $type=shift;
$self->{COORD_TYPE} = $type;
return $self;
}
##############################################################################
# set the coordinates "det" or "sky" to be used for WMAPs in spectra
# The default is "det".
##############################################################################
sub wmap_coord_type {
my $self=shift;
my $type=shift;
$self->{WMAP_COORD_TYPE} = $type;
return $self;
}
#############################################################################
# specify the output file name and type.
# Note that the mission and instrument must be specified before this
# method is called, since it needs to reference the MDB
#############################################################################
sub outfile {
my $self=shift;
my $file=shift;
my $type=shift || "event";
####################################
# clobber the outfile if it exists
# and remember its name
####################################
unlink $file;
$self->{OUTFILE}=$file;
####################################################
# set the parameters to set for this file type
####################################################
if( $type =~ /event/i ) {
#############################
# event files
#############################
$self->params({eventsout => $file,
imgfile => "NONE",
phafile => "NONE",
fitsbinlc => "NONE" });
} elsif($type =~ /image/i ) {
#######################################
# image
#######################################
$self->params({eventsout => "NONE",
imgfile => $file ,
phafile => "NONE",
fitsbinlc => "NONE" });
} elsif($type =~ /spec/i ) {
###########################################
# spectrum
###########################################
$self->params({eventsout => "NONE",
imgfile => "NONE",
phafile => $file ,
fitsbinlc => "NONE" });
} elsif($type =~ /light/i ) {
########################################
# light curve
#######################################
$self->params({eventsout => "NONE",
imgfile => "NONE",
phafile => "NONE",
fitsbinlc => $file });
} else {
##########################
# unknown
##########################
$self->setup_error(2,"Unknown file type $type for extractor");
}
#################################
# set generic parameters
#################################
$self->params({ecol => $self->mdb("ecol"),
ccol => $self->mdb("ccol"),
"gti" => $self->mdb("gti" ),
binf => $self->mdb("fbin"),
binh => $self->mdb("hbin"),
wtmapb => $self->mdb("wtmapb"),
wtmapfix => $self->mdb("wtmapfix"),
events => $self->mdb("events"),
copyall => "yes"});
####################################################
# get a sample fits input file which we can use to
# track down the column numbers for TLMIN/TLMAX
###################################################
unless($self->{SAMPLE_INPUT}) {
$self->log()->error(2,"Must Specify input files before output");
}
my $sample = Util::FITSfile->new($self->{SAMPLE_INPUT},
$self->mdb("events") );
####################################
# highest PHA channel
####################################
my $phamax = $self->mdb("phamax");
if($phamax eq "TLMAX") {
$phamax .= $sample->find_column($self->mdb("ecol"));
}
$self->params({phamax=>$phamax});
####################################
# filtering/image coordinates
####################################
if($self->{COORD_TYPE} =~ /det/i ) {
#######################
# detector coordinates
#######################
my $x_size = $self->mdb("detxsiz");
if($x_size eq "TLMAX") {
$x_size .= $sample->find_column($self->mdb("detx"));
}
my $y_size = $self->mdb("detysiz");
if($y_size eq "TLMAX") {
$y_size .= $sample->find_column($self->mdb("dety"));
}
$self->params({xcolf => $self->mdb("detx"),
ycolf => $self->mdb("dety"),
xfkey => $x_size,
yfkey => $y_size });
} elsif($self->{COORD_TYPE} =~ /raw/i ) {
#######################
# detector coordinates
#######################
my $x_size = $self->mdb("rawxsiz");
if($x_size eq "TLMAX") {
$x_size .= $sample->find_column($self->mdb("rawx"));
}
my $y_size = $self->mdb("rawysiz");
if($y_size eq "TLMAX") {
$y_size .= $sample->find_column($self->mdb("rawy"));
}
$self->params({xcolf => $self->mdb("rawx"),
ycolf => $self->mdb("rawy"),
xfkey => $x_size,
yfkey => $y_size });
} elsif($self->{COORD_TYPE} =~ /sky/i ) {
#######################
# sky coordinates
#######################
$self->params({xcolf => $self->mdb("x"),
ycolf => $self->mdb("y")});
}
####################################
# WMAP coordinates
####################################
if($self->{WMAP_COORD_TYPE} =~ /det/i ) {
#######################
# detector coordinates
#######################
$self->params({xcolh => $self->mdb("detx"),
ycolh => $self->mdb("dety") });
} elsif($self->{WMAP_COORD_TYPE} =~ /sky/i ) {
#######################
# detector coordinates
#######################
$self->params({xcolh => $self->mdb("x"),
ycolh => $self->mdb("y") });
}
return $self;
} # end of outfile method
###############
# METHODS:
###############
########################################################################
# simple error handler to report setup errors.
# It uses the log if available, but has a backup plan in case it isn't
#######################################################################
sub setup_error {
my $self=shift;
my $level=shift;
my $message=shift;
my $log=$self->log();
if($log) {
###########################
# use the log
###########################
$log->error($level,$message);
} else {
############################
# we got no log
############################
print STDERR "$message\n";
exit $level;
}
} # end of setup_error method
#######################################################################
# We override the Ftool run method because sometimes we will not want
# to run the tool at all - for example if time screening has excluded
# all times.
#######################################################################
sub run {
my $self = shift;
if($self->{NO_GOOD_TIME}) {
$self->log()->entry("No good time intervals. Extractor not run.");
return;
}
##############################################################
# if we get here then it is OK to actually run the extractor
#############################################################
$self->SUPER::run();
} # end of run method
######################################
# destructor - delete temporary files
######################################
sub DESTROY {
my $self=shift;
##############################
# delete the temporary files
##############################
unlink $self->{INFILES_TMP};
unlink $self->{REGION_TMP};
unlink $self->{GTI_TMP};
unlink $self->{ROW_SELECTION_FILE};
################################
# inherit the parent destructor
################################
$self->SUPER::DESTROY();
}