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