package ExpDetectPrep;
use strict;
use English;
use Carp;
use vars
qw(@ISA $VERSION $name $author $date $version @instList $numberOfStreams);
@ISA = qw(Module);
use ModuleResources;
use ModuleUtil;
use Regexp::Common qw(boolean);
use Astro::FITS::CFITSIO qw(TRUE FALSE);
$name = __PACKAGE__;
$VERSION = '0.26';
$version = $VERSION;
$author = 'Richard West,Dean Hinshaw,Duncan John Fyfe,Richard Saxton,Ian Stewart,Duncan Law-Green,Ed Chapin';
$date = '2013-09-05';
@instList = qw(emos1 emos2 epn);
my $xidBands = {'epn' => [qw(2 3 4)], 'emos1' => [qw(2 3 4)], 'emos2' => [qw(2 3 4)] };
sub numberOfStreams
{
1;
}
sub evaluateRules
{
my $inst = thisInstrument;
if ( $inst eq 'epn' )
{
return ignore()
if (
allIgnored(
module => 'MakePNImage'
, instrument => thisInstrument
)
);
start()
if (
allComplete(
module => 'MakePNImage'
, instrument => thisInstrument
)
);
}
else
{
return ignore()
if (
allIgnored(
module => 'MakeMOSImage'
, instrument => thisInstrument
)
);
start()
if (
allComplete(
module => 'MakeMOSImage'
, instrument => thisInstrument
)
);
}
}
sub performAction
{
info("Module version number: $version");
my @exp_id = listExposures(
instrument => thisInstrument
);
my @bandList = ( 1, 2, 3, 4, 5, 8 );
my %exposureLengths;
my (%filters, %modes);
foreach my $exp_id (@exp_id)
{
my @opt;
my $mode = getExposureProperty(
instrument => thisInstrument
, exp_id => $exp_id
, name => 'mode'
);
info("DEBUG: Exposure $exp_id Mode $mode");
if ( thisInstrument =~ /epn|mos/ )
{
my $filter = getExposureProperty(
instrument=> thisInstrument
, exp_id => $exp_id
, name => 'filter'
);
if ( $filter =~ /Open/ )
{
info("EPIC Open filter should not be used for source detection");
next;
}
}
if ( thisInstrument eq 'epn' )
{
my $evFile=findFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'EPIC PN imaging mode event list'
);
if ( $evFile && ModuleUtil::pnLowGainWarn($evFile) )
{
info("Exposure $exp_id: EPN exposure in low gain mode - won't perform source detection on it");
next;
}
}
elsif ( thisInstrument =~ /emos/ && $mode =~ /PrimePartialW6/i )
{
info("Exposure $exp_id: Will not perform source detection for MOS mode $mode");
next;
}
my $rawim = findFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, band => 1
, content => 'EPIC image'
, 'format' => 'FITS'
);
if ( !$rawim )
{
info("Exposure $exp_id: Can't find a band 1 image for this exposure - won't perform source detection on it");
next;
}
my $flareScreened;
@opt = ( file => $rawim , extension => 'PRIMARY' , keyword => 'BKGDSCRN');
if ( hasFITSKeyword( @opt ) )
{
$flareScreened = readFITSKeyword( @opt );
$flareScreened = ($flareScreened =~ /$RE{boolean}{true}/i)
? 1
: 0
;
}
if ( !$flareScreened )
{
info("Images were not flare screened - source detection not performed");
next;
}
my $filter = getExposureProperty(
instrument => thisInstrument
, exp_id => $exp_id
, name => 'filter'
);
unless ( $filter =~ /(Thin|Medium|Thick|Open)/ )
{
info("Could not match filter for energy conversion factor ($filter)");
next;
}
$filter = $1;
$filters{$exp_id} = $filter;
$modes{$exp_id} = $mode;
my $exposureLength = 0;
@opt = ( file => $rawim , extension => 'PRIMARY' , keyword => 'EXPOSURE' ) ;
if ( hasFITSKeyword( @opt ) )
{
$exposureLengths{$filter}{$mode} += readFITSKeyword( @opt );
}
}
return success() if (keys(%exposureLengths)<=0);
my $maxExposureLength = 0;
my ($chosenFilter, $chosenMode);
foreach my $filter (keys(%exposureLengths))
{
foreach my $mode (keys(%{$exposureLengths{$filter}}))
{
if ($exposureLengths{$filter}{$mode} > $maxExposureLength)
{
$maxExposureLength = $exposureLengths{$filter}{$mode};
$chosenFilter = $filter;
$chosenMode = $mode;
}
}
}
info("Chosen filter $chosenFilter");
info("Chosen mode $chosenMode");
my $tempset = newFile(
class => 'intermediate'
, content => 'Generic temporary dataset'
);
foreach my $band (@bandList)
{
my (@rawImList, @expImList,@expidList);
foreach my $exp_id (@exp_id)
{
next unless (($filters{$exp_id} eq $chosenFilter)
&& ( $modes{$exp_id} eq $chosenMode ) );
my $rawim = findFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, band => $band
, content => 'EPIC image'
, 'format' => 'FITS'
);
my $expim = findFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, band => $band
, content => 'EPIC exposure map'
, 'format' => 'FITS'
);
next unless ( $rawim && fileExists( file => $rawim )
&& $expim && fileExists( file => $expim )
);
push( @rawImList, $rawim );
push( @expImList, $expim );
push( @expidList, $exp_id );
}
unless ( @rawImList )
{
info("No images found");
next;
}
foreach my $exp_id (@expidList)
{
setExposureProperty( instrument => thisInstrument , exp_id => $exp_id , name => 'srcdet' , value => 0+TRUE );
}
info("Per camera merging of images for ",thisInstrument);
my $mergedRawim = newFile(
class => 'intermediate'
, instrument => thisInstrument
, band => $band
, content => 'merged image'
);
doCommand(
'imweightadd'
, imagesets => [@rawImList]
, withweights => 'no'
, outimageset => $mergedRawim
, tempset => $tempset
)
or return exception();
if (@rawImList > 1)
{
info("Calculating corrected ONTIME for merged images");
}
my $totalExposure;
my %ontime = ModuleUtil::ontimeFix( $mergedRawim , @rawImList );
if (my @ontime = sort( grep /ONTIME/,(keys %ontime)) )
{
my $num = scalar(@ontime);
my $numnum = $num+1;
my @attcomment = ('Sum of GTIs for a particular CCD') x $num;
push @attcomment , 'max of ONTIMEnn values' ;
my @wattcomment = ('Y') x $numnum;
my @atttype = ('real') x $numnum;
my @attname = (@ontime , 'EXPOSURE');
my @realval = (@ontime{@ontime} , $ontime{EXPOSURE});
doCommand('addattribute'
, set => $mergedRawim
,attributename => \@attname
,attributetype => \@atttype
,attributecomment => \@attcomment
,realvalue => \@realval
) or return exception();
doCommand('addattribute'
, set => $mergedRawim
,attributename => \@attname
,attributetype => \@atttype
,attributecomment => \@attcomment
,realvalue => \@realval
) or return exception();
$totalExposure = $ontime{EXPOSURE};
}
my @rawlistdtl = ( class => 'intermediate'
, content => 'raw image list'
, instrument => thisInstrument
, band => $band
, format => 'ASCII'
);
unless ( ModuleUtil::keepSafe ( file_identifier => \@rawlistdtl , text => \@rawImList ) )
{
info("Unable to store raw image list for later use.");
return exception();
}
foreach my $i (0 .. $#expidList)
{
my $exp_id = $expidList[$i];
my $kwd_str = 'EXPID'.sprintf("%2.2d", $i+1);
doCommand(
'addattribute', set => $mergedRawim.":PRIMARY"
, attributename => $kwd_str
, attributetype => 'string'
, stringvalue => $exp_id
)
or return exception();
}
my $mergedExpim = newFile(
class => 'intermediate'
, instrument => thisInstrument
, band => $band
, content => 'merged exposure map'
);
doCommand(
'imweightadd'
, imagesets => [@expImList]
, withweights => 'no'
, outimageset => $mergedExpim
, tempset => $tempset
) or return exception();
my @explistdtl = ( class => 'intermediate'
, content => 'exposure map list'
, instrument => thisInstrument
, band => $band
, format => 'ASCII'
);
unless ( ModuleUtil::keepSafe ( file_identifier => \@explistdtl , text => \@expImList ) )
{
info("Unable to store exposure map list for later use.");
return exception();
}
doCommand(
'addattribute', set => $mergedExpim.":PRIMARY"
, attributename => 'EXPOSURE'
, attributetype => 'real'
, realvalue => $totalExposure
) or return exception();
foreach my $i (0 .. $#expidList)
{
my $exp_id = $expidList[$i];
my $kwd_str = sprintf('EXPID%02d' , $i+1);
doCommand(
'addattribute', set => $mergedExpim.":PRIMARY"
, attributename => $kwd_str
, attributetype => 'string'
, stringvalue => $exp_id
) or return exception();
}
my (@flatExpImList, @imEventList);
foreach my $exp_id (@exp_id)
{
next unless ($filters{$exp_id} eq $chosenFilter
&& $modes{$exp_id} eq $chosenMode
);
my $flatExpIm = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, band => $band
, content => 'non-vig exposure map'
);
my $filteredList = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, band => $band
, content => 'image event list'
);
next unless $flatExpIm && $filteredList;
push( @flatExpImList, $flatExpIm );
push( @imEventList, $filteredList );
}
unless (@flatExpImList)
{
info("No flat expmaps found");
next;
}
my $mergedFlatExpim = newFile(
class => 'intermediate'
, instrument => thisInstrument
, band => $band
, content => 'merged non-vig exp map'
);
doCommand(
'imweightadd'
, imagesets => [@flatExpImList]
, withweights => 'no'
, outimageset => $mergedFlatExpim
, tempset => $tempset
)
or return exception();
doCommand(
'addattribute', set => $mergedFlatExpim.":PRIMARY"
, attributename => 'EXPOSURE'
, attributetype => 'real'
, realvalue => $totalExposure
) or return exception();
foreach my $i (0 .. $#expidList)
{
my $exp_id = $expidList[$i];
my $kwd_str = 'EXPID'.sprintf("%2.2d", $i+1);
doCommand(
'addattribute', set => $mergedFlatExpim.":PRIMARY"
, attributename => $kwd_str
, attributetype => 'string'
, stringvalue => $exp_id
) or return exception();
}
my $mergedFilteredList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, band => $band
, content => 'image merged event list'
);
my $imEventListSize = scalar(@imEventList);
info("DEBUG: + + + Number of image event lists: $imEventListSize");
my $fixtlm=1;
my ($tlmin6,$tlmax6,$tlmin7,$tlmax7);
if (scalar(@imEventList)>0) {
$tlmin6 = int(readFITSKeyword(file => $imEventList[0], extension => 'EVENTS', keyword => 'TLMIN6'));
$tlmax6 = int(readFITSKeyword(file => $imEventList[0], extension => 'EVENTS', keyword => 'TLMAX6'));
$tlmin7 = int(readFITSKeyword(file => $imEventList[0], extension => 'EVENTS', keyword => 'TLMIN7'));
$tlmax7 = int(readFITSKeyword(file => $imEventList[0], extension => 'EVENTS', keyword => 'TLMAX7'));
}
if (scalar(@imEventList)==1)
{
info("DEBUG: Branch 1");
$fixtlm=0;
copyFile(
source => $imEventList[0]
, destination => $mergedFilteredList
);
}
elsif(scalar(@imEventList)==2)
{
info("DEBUG: Branch 2");
if( ($tlmin6 != int(readFITSKeyword(file => $imEventList[1], extension => 'EVENTS', keyword => 'TLMIN6'))) ||
($tlmax6 != int(readFITSKeyword(file => $imEventList[1], extension => 'EVENTS', keyword => 'TLMAX6'))) ||
($tlmin7 != int(readFITSKeyword(file => $imEventList[1], extension => 'EVENTS', keyword => 'TLMIN7'))) ||
($tlmax7 != int(readFITSKeyword(file => $imEventList[1], extension => 'EVENTS', keyword => 'TLMAX7'))) )
{
$fixtlm=0;
info("Interesting: $imEventList[1] has different TL[MIN/MAX][6/7] keywords than $imEventList[0]. Won't fix merge output.");
}
doCommand(
'merge'
, set1 => $imEventList[0]
, set2 => $imEventList[1]
, outset => $mergedFilteredList
) or return exception();
} else {
info("DEBUG: Branch 3+");
my $prevFile=shift(@imEventList);
my $idx=0;
foreach my $file (@imEventList)
{
info("DEBUG: > > > Step $idx");
my $nextFile=newFile(class => 'intermediate',
content => 'Intermediate merged event list',
src_num => $idx
);
$idx++;
if( ($tlmin6 != int(readFITSKeyword(file => $file, extension => 'EVENTS', keyword => 'TLMIN6'))) ||
($tlmax6 != int(readFITSKeyword(file => $file, extension => 'EVENTS', keyword => 'TLMAX6'))) ||
($tlmin7 != int(readFITSKeyword(file => $file, extension => 'EVENTS', keyword => 'TLMIN7'))) ||
($tlmax7 != int(readFITSKeyword(file => $file, extension => 'EVENTS', keyword => 'TLMAX7'))) )
{
$fixtlm=0;
info("Interesting: $file has different TL[MIN/MAX][6/7] keywords than $imEventList[0]. Won't fix merge output.");
}
doCommand('merge',
,set1 => $prevFile
,set2 => $file
,outset => $nextFile
) or return exception();
$prevFile=$nextFile;
}
copyFile(source => $prevFile, destination => $mergedFilteredList);
}
if ($fixtlm) {
info("Merge may have broken TL[MIN/MAX}[6/7] keywords. Hacking $mergedFilteredList: TLMIN6=$tlmin6 TLMAX6=$tlmax6 TLMIN7=$tlmin7 TLMAX7=$tlmax7");
FITSUtils::writeKeyword($mergedFilteredList, "EVENTS", "TLMIN6" => {value => $tlmin6});
FITSUtils::writeKeyword($mergedFilteredList, "EVENTS", "TLMAX6" => {value => $tlmax6});
FITSUtils::writeKeyword($mergedFilteredList, "EVENTS", "TLMIN7" => {value => $tlmin7});
FITSUtils::writeKeyword($mergedFilteredList, "EVENTS", "TLMAX7" => {value => $tlmax7});
}
foreach my $i (0 .. $#expidList)
{
my $exp_id = $expidList[$i];
my $kwd_str = 'EXPID'.sprintf("%2.2d", $i+1);
doCommand(
'addattribute', set => $mergedFilteredList
, attributename => $kwd_str
, attributetype => 'string'
, stringvalue => $exp_id
) or return exception();
}
}
my @rawImList=();
my @expImList=();
my $inst = thisInstrument;
my @xidBandsThisInst = @{$xidBands->{$inst}};
foreach my $band (@xidBandsThisInst)
{
my $rawim = findFile(
class => 'intermediate'
, instrument => thisInstrument
, band => $band
, content => 'merged image'
);
my $expim = findFile(
class => 'intermediate'
, instrument => thisInstrument
, band => $band
, content => 'merged exposure map'
);
unless ( $rawim && fileExists( file => $rawim )
&& $expim && fileExists( file => $expim )
) {
info("Band $band image not found");
next;
}
push( @rawImList, $rawim );
push( @expImList, $expim );
}
unless (@rawImList == @xidBandsThisInst)
{
info("\nNo or missing XID-band constituent images.");
return success();
}
$tempset = newFile(
class => 'intermediate'
, content => 'Generic temporary dataset'
);
my $mergedRawim = newFile(
class => 'intermediate'
, instrument => thisInstrument
, band => 9
, content => 'merged image'
);
doCommand(
'imweightadd'
, imagesets => [@rawImList]
, withweights => 'no'
, outimageset => $mergedRawim
, tempset => $tempset
)
or return exception();
my $addedExpim = newFile(
class => 'intermediate'
, instrument => thisInstrument
, band => 9
, content => 'added exposure map'
);
my $mergedExpim = newFile(
class => 'intermediate'
, instrument => thisInstrument
, band => 9
, content => 'merged exposure map'
);
doCommand(
'imweightadd'
, imagesets => [@expImList]
, withweights => 'no'
, outimageset => $addedExpim
, tempset => $tempset
) or return exception();
my $div = @expImList;
doCommand('farith'
,infil1 => $addedExpim
,infil2 => $div
,outfil => $mergedExpim
,ops => 'DIV'
) or return exception();
return success();
}
1;