Download

# 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.04';
$version = $VERSION;
$author  = 'Ed Chapin, Jose Vicente Perea';
$date    = '2018-07-25';

#
# ChangeLog
# =========
#
# 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});

    }


    # Raise success flag
    return success();
}
1;