# After running pnEvents and MOSEvents, determine the TLMIN/TLMAX values that will encompass data from all
# EPIC cameras, and write these values to the imaging mode events files. This is initially aimed at allowing
# us to produce single images from all exposures within a mosaic.
package EPICImageSize;
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 POSIX; # for ceil
# Declare identity, version, author, date, etc.
$name = __PACKAGE__;
$VERSION = '0.09';
$version = $VERSION;
$author = 'Ed Chapin, Jose Vicente Perea';
$date = '2021-03-26';
#
# ChangeLog
# =========
#
# version 0.09 - 2021-03-26 JVP
# ------------
#
# + EPIC Merged Flare Background timeseries
# Final product moved from EPICEXP to EPICOBS group.
# Re-defined in Filemap as 'EPIC instrument flaring background'
#
# version 0.08 - 2021-03-23 JVP
# ------------
#
# + Merged flare bkg timeseries
# Correction: instrument definition in the plot filename
#
# version 0.07 - 2021-03-22 JVP
# ------------
#
# + Plot of the Merged flare bkg timeseries : X000FBKTSR*.PDF/PNG
# when there are more than 1 exposure per instrument
#
# version 0.06 - 2021-01-18 JVP
# ------------
#
# + Merged flare bkg timeseries to final product: X000FBKTSR
# when there are more than 1 exposure per instrument
#
# version 0.05 - 2020-12-03 JVP
# ------------
#
# + New module section:
# If more than one exposure per EPIC instrument, then
# Look for a common optimized flare background filter threshold from the merged timeseries
# of flare background timeseries from individual exposures.
# It has been included in this module because it runs once per all instruments/exposures
# just before making the EPIC images.
#
# version 0.04 - 2018-07-25 JVP
# ------------
#
# + Writing down the same Legal Limits ( tlmin6/7 , tlmax6/7 ) also into the FWC and OOT imaging mode event list
#
# version 0.03 - 2014-02-19 EC
# ------------
#
# + from PR: account for GTIs from preqgti, re-calculate tlmin/max from X,Y in events files
#
# version 0.02 - 2014-01-23 EC
# ------------
#
# + Remove ppd("000005") call when no imaging events found
#
# version 0.01 - 2013-02-06 EC
# ------------
#
# + First draft.
# Declare list of instruments this module is interested in
@instList = qw(epic);
# Number of streams
sub numberOfStreams
{
1;
}
#
sub evaluateRules
{
my $inst = thisInstrument;
# Ignore this module only if both the pnEvents and MOSEvents modules are completely ignored
return ignore()
if allIgnored( module => 'pnEvents' ) and allIgnored( module => 'MOSEvents' );
# Otherwise we should only start if all pnEvents and MOSEvents instances are either complete
# or ignored
start()
if (allComplete( module => 'pnEvents' ) and allComplete( module => 'MOSEvents' ));
}
#
sub performAction
{
info("Module version number: $version");
# Find all imaging event files
my @evFiles= ( findFile(class => 'product', content => 'EPIC PN imaging mode event list'),
findFile(class => 'product', content => 'EPIC MOS imaging mode event list') );
# Find attitude GTI file
my $attGTIFile = findFile(
class => 'intermediate'
, instrument => 'all'
, content => 'attitude GTI'
, required => 'true'
);
info("global Attitude GTI File: $attGTIFile");
# Return immediately if there are no event files
if ($#evFiles<0)
{
return success();
}
# load in the GTI data
my $gti = readFITSTable( file => $attGTIFile,
extension => "STDGTI",
colname => [qw( START STOP )] );
# find absolute minimum and maximum data ranges of X,Y columns,
# TLMIN/TLMAX, defining the actual image size
my ($tmin6, $tmax6, $tmin7, $tmax7);
foreach my $evFile (@evFiles)
{
my $ev = readFITSTable( file => $evFile,
extension => "EVENTS",
colname => [qw( TIME X Y )] );
my $xnull = readFITSKeyword( file => $evFile
,extension => 'EVENTS'
,keyword => 'TNULL6');
my $ynull = readFITSKeyword( file => $evFile
,extension => 'EVENTS'
,keyword => 'TNULL7');
info("X,Y NULL values: $xnull , $ynull\n");
my $xmin;
my $xmax;
my $ymin;
my $ymax;
my $init=1; # unset once values are set to something
# Loop over events
for my $j ( 0 .. $#{$ev->{TIME}} ) {
my $time = ${$ev->{TIME}}[$j];
my $x = ${$ev->{X}}[$j];
my $y = ${$ev->{Y}}[$j];
# Loop over GTIs
for my $i ( 0 .. $#{$gti->{START}} ) {
# Is the event within the GTI?
if ( $time >= ${$gti->{START}}[$i] &&
$time <= ${$gti->{STOP}}[$i] &&
$x != $xnull &&
$y != $ynull ) {
if ( $init ) {
$xmin = $x;
$xmax = $x;
$ymin = $y;
$ymax = $y;
undef $init;
info("Initial: $xmin $xmax $ymin $ymax");
} else {
if ($x < $xmin) { $xmin = $x; }
if ($x > $xmax) { $xmax = $x; }
if ($y < $ymin) { $ymin = $y; }
if ($y > $ymax) { $ymax = $y; }
unless( defined($xmin) && defined($xmax) && defined($ymin) && defined($ymax) ) {
info("Huh? Undefined bounds: xmin=$xmin xmax=$xmax ymin=$ymin ymax=$ymax");
}
}
}
}
}
info("First pass $evFile: xmin=$xmin xmax=$xmax ymin=$ymin ymax=$ymax");
# Calculate quantized values
my $thistmin6;
my $thistmax6;
my $thistmin7;
my $thistmax7;
if( defined($xmin) && defined($xmax) && defined($ymin) && defined($ymax) ) {
my $xminscaled = ceil(-1.0*$xmin/3600.0);
my $xmaxscaled = ceil(($xmax-51840)/3600.0);
my $yminscaled = ceil(-1.0*$ymin/3600.0);
my $ymaxscaled = ceil(($ymax-51840)/3600.0);
$thistmin6= 1-3600*( ($xminscaled > 0) ? $xminscaled : 0 );
$thistmax6 = 51840+3600*( ($xmaxscaled > 0) ? $xmaxscaled : 0);
$thistmin7= 1-3600*( ($yminscaled > 0) ? $yminscaled : 0 );
$thistmax7 = 51840+3600*( ($ymaxscaled > 0) ? $ymaxscaled : 0);
# the previous code should be doing the following pseudocode:
#$thistmin6= 1-3600*max(0,ceil(-1.0*$xmin/3600.0)) ;
#$thistmax6 = 51840+3600*max(0,ceil(($xmax-51840)/3600.0));
#$thistmin7 = 1-3600*max(0,ceil(-1.0*$ymin/3600.0));
#$thistmax7 = 51840+3600*max(0,ceil(($ymax-51840)/3600.0));
} else {
info("No valid events during GTIs found in $evFile.");
}
# Old method based on FITS keywords
#my $thistmin6 = int(readFITSKeyword(file => $evFile, extension => 'EVENTS', keyword => 'TLMIN6'));
#my $thistmax6 = int(readFITSKeyword(file => $evFile, extension => 'EVENTS', keyword => 'TLMAX6'));
#my $thistmin7 = int(readFITSKeyword(file => $evFile, extension => 'EVENTS', keyword => 'TLMIN7'));
#my $thistmax7 = int(readFITSKeyword(file => $evFile, extension => 'EVENTS', keyword => 'TLMAX7'));
info("First pass $evFile: tmin6=$thistmin6 tmax6=$thistmax6 tmin7=$thistmin7 tmax7=$thistmax7");
if ($evFile eq $evFiles[0])
{
# initialize values if this is the first file
$tmin6 = $thistmin6;
$tmax6 = $thistmax6;
$tmin7 = $thistmin7;
$tmax7 = $thistmax7;
}
else
{
# otherwise check to see if we need to extend the extrema
$tmin6 = $thistmin6 if $thistmin6 < $tmin6;
$tmax6 = $thistmax6 if $thistmax6 > $tmax6;
$tmin7 = $thistmin7 if $thistmin7 < $tmin7;
$tmax7 = $thistmax7 if $thistmax7 > $tmax7;
}
info("Current max ranges: $tmin6 $tmax6 $tmin7 $tmax7");
}
unless( defined($tmin6) && defined($tmax6) && defined($tmin7) && defined($tmax7) ) {
info("No valid events found, so TLMIN/MAX will not be updated");
return success();
}
info("Full set of EPIC events encompassed by TLMIN6=$tmin6 TLMAX6=$tmax6 TLMIN7=$tmin7 TLMAX7=$tmax7");
# Loop again over events files and write out the extra for the maximum legal ranges TLMIN/TLMAX
foreach my $evFile (@evFiles)
{
# addattribute won't let us write TLMIN/MAX keywords as they are special keywords
#doCommand( 'addattribute'
# ,attributename => [qw(TLMIN6 TLMAX6 TLMIN7 TLMAX7)]
# ,attributetype => [qw(integer integer integer integer)]
# ,integervalue => [ $tmin6, $tmax6, $tmin7, $tmax7 ]
# , set => ["$evFile:EVENTS"]
# ) or exception();
info("Writing TLMIN/MAX{6,7} to $evFile. TLMIN6=$tmin6 TLMAX6=$tmax6 TLMIN7=$tmin7 TLMAX7=$tmax7");
FITSUtils::writeKeyword($evFile, "EVENTS", "TLMIN6" => {value => $tmin6});
FITSUtils::writeKeyword($evFile, "EVENTS", "TLMAX6" => {value => $tmax6});
FITSUtils::writeKeyword($evFile, "EVENTS", "TLMIN7" => {value => $tmin7});
FITSUtils::writeKeyword($evFile, "EVENTS", "TLMAX7" => {value => $tmax7});
}
# Loop also over imaging events files for FWC and OOT
my @evFilesFWCandOOT= ( findFile(class => 'intermediate', content => 'EPIC PN imaging mode FWC event list'),
findFile(class => 'intermediate', content => 'EPIC PN imaging mode OOT event list'),
findFile(class => 'intermediate', content => 'EPIC MOS imaging mode FWC event list') );
foreach my $evFile (@evFilesFWCandOOT)
{
info("Writing TLMIN/MAX{6,7} to $evFile. TLMIN6=$tmin6 TLMAX6=$tmax6 TLMIN7=$tmin7 TLMAX7=$tmax7");
FITSUtils::writeKeyword($evFile, "EVENTS", "TLMIN6" => {value => $tmin6});
FITSUtils::writeKeyword($evFile, "EVENTS", "TLMAX6" => {value => $tmax6});
FITSUtils::writeKeyword($evFile, "EVENTS", "TLMIN7" => {value => $tmin7});
FITSUtils::writeKeyword($evFile, "EVENTS", "TLMAX7" => {value => $tmax7});
}
###########################################################################################
# Looking for a common flare background filtering threshold for cases with #exposures > 1
#
foreach my $inst (qw(epn emos1 emos2)){
my @FBKTSRFiles = findFile(
class => 'product'
, instrument => $inst
#, exp_id => $exp_id
, content => 'EPIC flare background timeseries'
, format => 'FITS'
);
foreach my $f (@FBKTSRFiles){
if( ! hasFITSKeyword(file => $f, extension => 'RATE', keyword => 'BACKSCAL') ){
info("Wrong EPIC flare background timeseries $f. BACKSCAL keyword not found. File excluded from timeseries list.");
@FBKTSRFiles = grep {$_ ne $f} @FBKTSRFiles;} # File $f removed from @FBKTSRFiles list
}
if ( scalar(@FBKTSRFiles) > 1){
info("There are more than one exposure for $inst.");
info("Looking for a common optimized flare background filter threshold from the merged timeseries");
my $mergedBackRates = newFile(
class => 'product'
, instrument => $inst
, content => 'EPIC instrument flaring background'
, format => 'FITS'
);
# No exp_id set => merged timeseries (X000)
my $ratesFiles = join(',', @FBKTSRFiles);
info("Merging timeseries: $ratesFiles");
doCommand(
'mergeEPICflareBkgTS.py'
, ratesFiles => $ratesFiles
, mergedRates => $mergedBackRates
) or return exception();
info("Merged timeseries: $mergedBackRates");
my $snDiagFile = newFile(
class => 'intermediate'
, instrument => $inst
, content => 'Merged Flare GTI SN diagnostic'
, band => 8
, format => 'FITS'
);
my $retval = 0;
my $rcut = `bkgoptrate -V 0 -w 0 tssettabname=${mergedBackRates}:RATE doinsertkwds=yes dooutputsntab=yes ignorelowcnttail=yes snsettabname=${snDiagFile}:SN_TO_BKGCUT`;
info("CMD: bkgoptrate -V 0 -w 0 tssettabname=${mergedBackRates}:RATE doinsertkwds=yes dooutputsntab=yes ignorelowcnttail=yes snsettabname=${snDiagFile}:SN_TO_BKGCUT");
$retval = $?;
if ($retval ne 0) { return exception("Error $retval returned by bkgoptrate: $!"); }
my $optFlareGti;
if ($rcut eq '') {
info("bkgoptrate call returns no value, skipping optimised flare GTI generation for merged timeseries...");
} else {
info("bkgoptrate returns optimum cut threshold: $rcut");
# append SN_TO_BKGCUT table to FBKTSR product
doCommand(
'ftappend'
, infile => "${snDiagFile}[SN_TO_BKGCUT]"
, outfile => $mergedBackRates
) or return exception();
# add FLCUTTHR header value
doCommand(
'addattribute'
, set => "${mergedBackRates}:RATE"
, attributename => 'FLCUTTHR'
, attributetype => 'real'
, attributecomment => '"\"Optimised flare cut threshold\""'
, realvalue => $rcut
) or return exception();
# Now generate the optimised flare GTI files
$optFlareGti = newFile(
class => 'intermediate'
, instrument => $inst
, content => "Optimised flare GTI from merged timeseries"
, band => 8
, format => 'FITS'
);
info("Creating GTI files with this threshold: $rcut");
doCommand(
'tabgtigen'
, table => $mergedBackRates
, gtiset => $optFlareGti
, expression => "RATE < ${rcut}"
, mingtisize => 100
) or return exception();
info("Make full in-band, source-free merged light curve plots");
#
if ( hasFITSExtension(file => $mergedBackRates, extension => "RATE") &&
hasFITSExtension(file => $mergedBackRates, extension => "SN_TO_BKGCUT") ){
my $flareBkgPDF0 = newFile(
class => 'product'
, instrument => $inst
, content => 'EPIC instrument flaring background'
, format => 'PDF'
);
my $flareBkgPNG0 = newFile(
class => 'product'
, instrument => $inst
, content => 'EPIC instrument flaring background'
, format => 'PNG'
);
doCommand(
'plotFlrBkgTimeSeries.py'
, rateset => $mergedBackRates
, plotfile => $flareBkgPNG0
, pdfplotfile => $flareBkgPDF0
, plotformat => "both"
) or return exception();
} else {
info("No RATE or SN_TO_BKGCUT extensions found in source-free merged light curve. No plot created.");
}
}
}
}
# Raise success flag
return success();
}
1;