Download

# 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;