# Library of utility FUNCTIONS common to several modules.
package ModuleUtil;
use strict;
use English;
use Carp;
use ModuleResources;
use FITSUtils;
use File::Temp;
use Regexp::Common;
use vars qw( $VERSION );
$VERSION = 1.0;
#----------------------------
# output pattern of fimgstat
#----------------------------
#The sum of the selected image = 118590.000000
#The mean of the selected image = 0.282422
#The standard deviation of the selected image = 19.612485
#The number of points used in calculation = 419904
#The minimum of selected image = 0.0
#The maximum of selected image = 5573.0
#The location of minimum is at pixel number = (1,1)
#The location of maximum is at pixel number = (325,322)
use vars qw( @Fimgstat_Patterns );
# Using an array rather than a hash because the input data is ordered
# so this will be quicker.
@Fimgstat_Patterns = (
sum => qr/^The sum of the selected image\s*?=\s*?$RE{num}{real}{-keep}\s*$/
,mean => qr/^The mean of the selected image\s*?=\s*?$RE{num}{real}{-keep}\s*$/
,stddev => qr/^The standard deviation of the selected image\s*?=\s*?$RE{num}{real}{-keep}\s*$/
,count => qr/^The number of points used in calculation\s*?=\s*?$RE{num}{int}{-keep}\s*?/
,min => qr/^The minimum of selected image\s*?=\s*?$RE{num}{real}{-keep}\s*$/
,max => qr/^The maximum of selected image\s*?=\s*?$RE{num}{real}{-keep}\s*$/
,minpos => qr/^The location of minimum is at pixel number\s*?=\s*?\(($RE{num}{int},$RE{num}{int})\)\s*$/
,maxpos => qr/^The location of maximum is at pixel number\s*?=\s*?\(($RE{num}{int},$RE{num}{int})\)\s*$/
);
#------------------------------
# output pattern of fstatistic
#------------------------------
# The sum of the selected column is 0.0000000
# The mean of the selected column is 0.0000000
# The standard deviation of the selected column is 0.0000000
# The minimum of selected column is 0.0000000
# The maximum of selected column is 0.0000000
# The number of points used in calculation is 4347
use vars qw( @Fstatistic_Patterns );
# Using an array rather than a hash because the input data is ordered
# so this will be quicker.
@Fstatistic_Patterns = (
sum => qr/^\s*The sum of the selected column is\s*$RE{num}{real}{-keep}\s*$/
,mean => qr/^\s*The mean of the selected column is\s*$RE{num}{real}{-keep}\s*$/
,stddev => qr/^\s*The standard deviation of the selected column is\s*$RE{num}{real}{-keep}\s*$/
,min => qr/^\s*The minimum of selected column is\s*$RE{num}{real}{-keep}\s*$/
,max => qr/^\s*The maximum of selected column is\s*$RE{num}{real}{-keep}\s*$/
,count => qr/^\s*The number of points used in calculation is\s*?$RE{num}{int}{-keep}\s*?/
);
sub isSourceOnChip
{
# Use echeckregion to determine if a source is within the detector window or not.
# Sometimes a source extraction region clips a chip giving enough events
# to pass the .non-zero' events test but there is no RMF as the source is
# in a void.
# Returns : 1 = onchip, 0 = offchip
my ($srcexp , $table , $output) = @_;
my $status = &doCommand(
'echeckregion'
, srcexp => $srcexp
, table => $table
, output => $output
) or exception();
my $rtn = 0;
if ( &fileExists( file => $output ) )
{
my $fh;
open($fh , "<$output")
or Exception->error('FILESYSTEM','Unable to open echeckregion output file.' , $output )
;
my @tmp = grep /Status\s+of\s+the\s+centre\s+of\s+the\s+input\s+region\s+is:/,(<$fh>);
close($fh);
if (@tmp)
{
$tmp[0] =~ /Status\s+of\s+the\s+centre\s+of\s+the\s+input\s+region\s+is:\s*(\d+)/;
$rtn = $1;
}
# else assume the answr is no (0).
}
else
{
Msg->info("Unable to find output from echeckregion.");
exception();
}
return $rtn;
}
sub parseFimgstatOutput
{
my $f = shift;
Exception->error('BADMETHODCALL','Called with non-existant file: ' , $f )
unless ( -f $f )
;
# Search the output from fimgstat for what we want
my $fh;
open($fh , "<$f")
or Exception->error('FILESYSTEM','Unable to open fimgstat output file.' , $f )
;
my @tmpstats = <$fh>;
close($fh);
my %h;
my $i = 0;
foreach my $line (@tmpstats)
{
for ( ; $i< @Fimgstat_Patterns ; $i+=2 )
{
my @q = $line =~ $Fimgstat_Patterns[$i+1];
next unless @q;
if ($q[0] !~ /,/)
{
$h{ $Fimgstat_Patterns[$i] } = $q[0];
}
else
{
$h{ $Fimgstat_Patterns[$i] } = [ split /,/,$q[0] ];
}
last;
}
}
return %h;
}
sub getImageStats
{
# Call fimgstat and put the results into a hash
my ($imgFile , $statFile) = @_;
$statFile ||= File::Temp->new()->filename();
unless ( -f $statFile)
{
&doCommand (
'fimgstat', infile => "$imgFile"
, threshlo => 'INDEF'
, threshup => 'INDEF'
, outfile => "$statFile"
, clobber => 'Y'
, mode => 'ql'
, sum => 0.0
, mean => 0.0
, rms => 0.0
, num => 0
, min => 0.0
, max => 0.0
, xmin => 0
, ymin => 0
, xmax => 0
, ymax => 0
) or exception()
;
}
my %h = &parseFimgstatOutput( $statFile );
return %h;
}
sub isImageEmpty
{
# get context for error reporting.
#
# Use fimgstat to determine if a file has sufficient range for implot.
# Return 1 = empty |0 = not empty
#
my ($imgFile , $statFile) = @_;
my $isEmpty = 0;
my %h = &getImageStats($imgFile,$statFile);
unless (defined $h{max} && defined $h{min} )
{
# We didn't find what we expected to find.
Exception->error('FIMGSTAT','Unable to parse fimgstat output.');
}
if ( !$h{max} || ( $h{max} - $h{min} ) <= 1.0 )
{
Msg->info("The image $imgFile is empty.");
$isEmpty = 1;
}
return $isEmpty;
}
sub isImageLimit
{
#
# Use fimgstat to determine if image has sufficient counts for subsequent analysis
# Return 1 = below limit |0 = equal or above limit
# $imgLimit is number of limiting counts
my ($imgFile , $statFile, $imgLimit) = @_;
my $isLimit = 0;
my %h = &getImageStats($imgFile,$statFile);
unless (defined $h{max} && defined $h{min} )
{
# We didn't find what we expected to find.
Exception->error('FIMGSTAT','Unable to parse fimgstat output.');
}
if ( $h{sum} < $imgLimit )
{
my $sum = $h{sum};
Msg->info("The image $imgFile has $sum total counts, below limit $imgLimit");
$isLimit = 1;
}
return $isLimit;
}
sub insertdb
{
# Attatch an arbitrary data structure to the sequence database.
# The structure can be retrieved using the given key.
# Once created given key,value pair become read only.
# %d = ()
my ($key,$data) = (@_);
# Turn an array reference key into a usable key
if (UNIVERSAL::isa($key,'ARRAY') )
{
$key = join(':',@$key);
}
if (exists $ModuleResources::Seqdb->{_moduledb}{$key})
{
Exception->error("Attempt to amend read-only module database entry");
}
$ModuleResources::Seqdb->{_moduledb}{$key} = Storable::dclone($data);
return $key;
}
sub selectdb
{
# Attatch an arbitrary data structure to the sequence database.
# The structure can be retrieved using the given key.
# Once created given key,value pair become read only.
# %d = ()
my ($key,$data) = (@_);
# Turn an array reference key into a usable key
if ( UNIVERSAL::isa($key,'ARRAY') )
{
$key = join(':',@$key);
}
if (exists $ModuleResources::Seqdb->{_moduledb}{$key})
{
return $ModuleResources::Seqdb->{_moduledb}{$key};
}
return undef;
}
sub keepSafe
{
# %d = ( file_identifier => [] , text => [] )
# fileid == enough to find an existing file or create a new one.
my %d = (@_);
my $rtn = 0;
my $safe = findFile( @{$d{file_identifier}} );
unless ( $safe && fileExists($safe) )
{
$safe = newFile( @{$d{file_identifier}} );
}
if ($safe)
{
my @d = @{ $d{text} };
foreach my $d (@d) { $d .= "\n"; }
writeASCIIFile( name => $safe
, text => \@d
, append => 1
);
$rtn = 1;
}
return $rtn;
}
sub ontimeFix
{
my ($output , @image) = (@_);
my @ontime = qw(
ONTIME01
ONTIME02
ONTIME03
ONTIME04
ONTIME05
ONTIME06
ONTIME07
ONTIME08
ONTIME09
ONTIME10
ONTIME11
ONTIME12
);
Msg->mark();
my %ontime = ();
my $exposure = 0;
foreach my $img (@image)
{
next unless -f $img;
foreach my $kw ( @ontime )
{
my @opt = ($img,'PRIMARY',$kw);
next unless FITSUtils::hasKeyword(@opt);
my $v = FITSUtils::readKeyword(@opt);
$ontime{$kw} += FITSUtils::readKeyword(@opt);;
if ($exposure < $ontime{$kw})
{
$exposure = $ontime{$kw};
}
}
}
unless (keys %ontime)
{
Msg->info("No ONTIME values in input images.");
}
return (%ontime , EXPOSURE => $exposure);
}
sub pnLowGainWarn
{
my ( $file , $h ) = @_;
my @ext = qw(
EXPOSU01
EXPOSU02
EXPOSU03
EXPOSU04
EXPOSU05
EXPOSU06
EXPOSU07
EXPOSU08
EXPOSU09
EXPOSU00
EXPOSU11
EXPOSU12
);
my @warn;
foreach my $ext ( @ext )
{
next unless ( FITSUtils::hasExtension( $file , $ext ) );
foreach my $kw (qw( GAINCAME WARN_106 ))
{
my @opt = ( $file , $ext , $kw );
if ( FITSUtils::hasKeyword( @opt ) )
{
$h->{$ext}{$kw} = FITSUtils::readKeyword(@opt);
}
}
push @warn , $ext if ( ($h->{$ext}{GAINCAME} && $h->{$ext}{GAINCAME} eq 'LOW' )
&& ( $h->{$ext}{WARN_106} && $h->{$ext}{WARN_106} =~ /CCDlowGain/i )
);
}
return @warn;
}
sub select_rgssources
{
my ($file , $extension , $haveEpicList , $rate_cutoff ) = (@_);
# Select additional EPIC sources for spectral extraction from the RGS source list.
# Returns a source selection expression for passing to rgsregions.
# Also works
my $primesrc = readFITSKeyword(
file => $file
, extension => $extension
, keyword => "PRIMESRC"
);
my @idx;
# Only if we are playing with an EPIC source list.
# Doesn't make sense otherwise.
if ( $haveEpicList )
{
my $data = readFITSTable(
file => $file
, extension => $extension
, colname => [qw( INDEX LABEL RATE FLAG )]
);
my $processables = 0;
my $hasepic = 0;
for my $i ( 0 .. $data->{-nrows} -1 )
{
# Skipables.
# avoid where:
# FLAG set
# avoid duplicating the PRIMESRC
# too faint.
# non-EPIC source (PROPOSAL, ONAXIS , WHOLEFIELD etc)
# Make sure it is processable.
next if ( ( $data->{FLAG}[$i] ne 0 )
|| ( $data->{LABEL}[$i] !~ /EPIC/ )
);
# Record the fact we have any EPIC sources at all.
$hasepic ||= 1;
next unless ( $data->{RATE}[$i] > $rate_cutoff );
# Count EPIC sources which pass the above.
$processables++;
# but avoid adding EPIC PRIMESRC to the list.
next if ( $data->{INDEX}[$i] eq $primesrc );
push @idx, $i;
}
# Sort the list of indexes on RATE
@idx = sort { $data->{RATE}[$a] <=> $data->{RATE}[$b] } @idx;
#convert @idx to INDEX column values.
@idx = map { $data->{INDEX}[$_]; } @idx;
doCommand('addattribute'
, set => "$file:SRCLIST"
, attributename => [ 'EPMINRTE' , 'EPRGSRCN' ]
, attributetype => [ 'real' , 'integer' ]
, attributecomment => [
'EPIC source RATE Threshold'
,'Number of EPIC source within RGS FOV and RATE > EPMINRTE and FLAG 0'
]
,integervalue => $processables
,realvalue => $rate_cutoff
) or exception();
}
my @rtn = (@idx)
? ( "INDEX==#PRIMESRC || INDEX==$idx[-1]" , $primesrc , $idx[-1] )
: ( "INDEX==#PRIMESRC" , $primesrc )
;
return @rtn;
}
sub get_rgssources_param
{
my ($file, $extension , $instrument) = (@_);
# Select sources from RGS source list so details can be recorded in the source spectra
my $primesrc = readFITSKeyword( file => $file
, extension => $extension
, keyword => "PRIMESRC"
);
my %rtn;
my $data = readFITSTable(
file => $file
, extension => $extension
, colname => [qw( INDEX LABEL RA DEC RATE PROCESS FLAG )]
);
my $processables = 0;
my $hasepic = 0;
for my $i ( 0 .. $data->{-nrows} -1 )
{
# Skipables.
# avoid where:
# FLAG set
# avoid duplicating the PRIMESRC
# too faint.
# non-EPIC source (PROPOSAL, ONAXIS , WHOLEFIELD etc)
# Make sure it is processable.
next if ( ( $data->{PROCESS}[$i] =~ /$RE{boolean}{false}/i )
|| ( ! hasFITSExtension( file => $file
, extension => "${instrument}_SRC".$data->{INDEX}[$i]."_SPATIAL"
)));
# Keep EPIC source details to be recorded in the spectrum headers.
$rtn{ $data->{INDEX}[$i] } = {
INDEX => $data->{INDEX}[$i]
,LABEL => $data->{LABEL}[$i]
,RA => $data->{RA}[$i]
,DEC => $data->{DEC}[$i]
,RATE => $data->{RATE}[$i]
,PROCESS => $data->{PROCESS}[$i]
,FLAG => $data->{FLAG}[$i]
};
}
return %rtn;
}
sub parseFstatisticOutput
{
my $f = shift;
Exception->error('BADMETHODCALL','Called with non-existant file: ' , $f )
unless ( -f $f )
;
# Search the output from fstatistic for what we want
my $fh;
open($fh , "<$f")
or Exception->error('FILESYSTEM','Unable to open fstatistic output file.' , $f )
;
my @tmpstats = <$fh>;
close($fh);
my %h;
my $i = 0;
foreach my $line (@tmpstats)
{
for ( ; $i< @Fstatistic_Patterns ; $i+=2 )
{
my @q = $line =~ $Fstatistic_Patterns[$i+1];
next unless @q;
if ($q[0] !~ /,/)
{
$h{ $Fstatistic_Patterns[$i] } = $q[0];
}
else
{
$h{ $Fstatistic_Patterns[$i] } = [ split /,/,$q[0] ];
}
last;
}
}
return %h;
}
sub isColumnEmpty
{
# get context for error reporting.
#
# Use fstatistic to determine if a table column is identically 0
# Return 1 = empty |0 = not empty
#
my ($inFile , $colName, $statFile) = @_;
my $isEmpty = 0;
$statFile ||= File::Temp->new()->filename();
unless ( -f $statFile)
{
&doCommand (
'fstatistic', infile => "$inFile"
, colname => "$colName"
, rows => "-"
, outfile => "$statFile"
) or exception()
;
}
my %h = &parseFstatisticOutput( $statFile );
unless (defined $h{max} && defined $h{min} )
{
# We didn't find what we expected to find.
Exception->error('FSTATISTIC','Unable to parse fstatistic output.');
}
if ( ($h{max}==0) && ($h{min}==0) )
{
Msg->info("The column $colName of $inFile is identically 0.");
$isEmpty = 1;
}
return $isEmpty;
}
# Modules like ExpDetect were not originally set up for mosaics,
# although in principle they might be OK if the individual pointings
# are all very close to one another. This subroutine simply checks a
# reference image to determine whether its angular extent is larger
# than 0.82 deg along either axis as a simple test of whether it is a
# large mosaic or not.
sub isSmallMap
{
my $imFile = $_[0];
my $status = 0;
return $status unless $imFile;
#my $cunit1 = readFITSKeyword(file => $imFile, extension => 1, keyword => 'CUNIT1');
#my $cunit2 = readFITSKeyword(file => $imFile, extension => 1, keyword => 'CUNIT2');
#
#Exception->error('Module', 'Supplied image file has the wrong axis units:', $imFile)
# unless ( ($cunit1 =~ m/deg/i) and ($cunit2 =~ m/deg/i) );
info( "Warning: not checking units since some maps are lacking CUNIT keywords (e.g., expposure maps)");
my $xsize = readFITSKeyword(file => $imFile, extension => 1, keyword => 'CDELT1') *
readFITSKeyword(file => $imFile, extension => 1, keyword => 'NAXIS1');
my $ysize = readFITSKeyword(file => $imFile, extension => 1, keyword => 'CDELT2') *
readFITSKeyword(file => $imFile, extension => 1, keyword => 'NAXIS2');
$xsize = abs($xsize);
$ysize = abs($ysize);
$status = 1 if ($xsize < 0.82) and ($ysize < 0.82);
info("Dimensions of the test image $imFile are $xsize x $ysize. isSmallMap=$status.");
return $status;
}
1;