package pnEvents;
use strict;
use English;
use Carp;
use List::Util qw(max min);
use Regexp::Common qw(boolean);
use vars
  qw(@ISA $VERSION $name $author $date $version @instList $numberOfStreams);
@ISA = qw(Module);
use ModuleResources;

# Declare identity, version, author, date, etc.
$name    = __PACKAGE__;
$VERSION = '2.72';
$version = $VERSION;
$author  = 'Dean Hinshaw,Richard West,Duncan John Fyfe,Ian Stewart,Duncan Law-Green,Eduardo Ojero,Ed Chapin, Jose Vicente Perea';
$date    = '2021-03-22';

#
# ChangeLog
# =========
#
# Version 2.72 - 2021-03-22 (JVP)
# ------------
#
# + EPIC Flare Background timeseries plot: source free background timeseries and Signal-to-noise plots
#   Old code removed (fplot, gs and psnup) and replaced by Python code : plotFlrBkgTimeSeries.py
#   (Besides PDF is also created as PNG file)
#
# Version 2.71 - 2021-01-18 (JVP)
# ------------
#
# + Change back to "pswrite" from "ps2write"
#   It must be changed again in future OS with newer version of ghostscript
#
# Version 2.70 - 2020-11-27 (JVP)
# ------------
#
# + Ghostscript device deprecated: "pswrite", changed to "ps2write"
#   (Option removed in ghostscript-9.10)
#
# Version 2.69 - 2019-04-10 (JVP)
# ------------
#
# + Skip SFBLC exposure maps if SFBLC input image has EXPOSURE = 0
#   'eexpmap' task crashes otherwise
#
# Version 2.68 - 2018-07-25 (JVP)
# ------------
#
# + Writing down the same Legal Limits ( tlmin6/7 , tlmax6/7 ) also in the FWC imaging event list: $finFWCEvents
#   Be aware that TLMIN6, TLMAX6, TLMIN7, TLMAX7 are again re-written into the final products,
#   FWC and OOT imaging event files by the module 'EPICImageSize.pm'
#
# Version 2.67 - 2018-07-09 (JVP)
# ------------
#
# + HKAUX added to othertables option on 'evlistcomb' call
# + Default param values for backgroundtres=Y and backgroundtbin=100 explicitly set on 'epevents' call
#
# Version 2.66 - 2018-06-14 (JVP)
# ------------
#
# + Small comment
#
# Version 2.65 - 2018-06-05 (JVP)
# ------------
#
# + Search for FWC data (evqpb) only if the observing mode is 'PrimeFullWindow' or 'PrimeFullWindowExtended'
#
# Version 2.64 - 2017-11-17 (JVP)
# ------------
#
# + Add epexposure call w/time randomisation for OoT events only in Imaging mode
#
# Version 2.63 - 2017-09-21 (JVP)
# ------------
#
# + Include randomizatoin also for '$tmpootEvents' (epexposure)
# + Writing down the same Legal Limits ( tlmin6/7 , tlmax6/7 ) in both $finEvents and $finootEvents
#   to ensure the same exactly borders for both event files.
# + Note about $attFile
#     Same $attFile is applied for $evFile2/$$evootFile2, $finEvents/$finootEvents and $attGTIFilteredEvents/attGTIFilteredootEvents
#     to ensure the same astrometry (sky projection, etc) in both event files: all and oot events.
#
# Version 2.62 - 2017-09-19 (JVP)
# ------------
#
# + Out-of-Time and FWC correction (first steps in PN)
#
# Version 2.61 - 2017-01-24 (JVP)
# ------------
#
# + The Mask is produced although image counts < 1000.
#   SFBLCBackScal calculated although image counts < 1000. Calculated only from the Mask.
#
# Version 2.60 - 2016-12-16 (JVP)
# ------------
#
# + BACKSCAL (AREASCALE parameter) included as keyword in the 'EPIC FLARE BACKGROUND TIMESERIES'
#   Image pixel size defined at the beginning as x/ybin_SFBLC = 80
#
# Version 2.59 - 2016-12-14 (JVP)
# ------------
#
# + Calculate AREASCALE parameter in background estimate for flaring filtering
#   New SAS task included: regionmask
#
# Version 2.58 - 2016-03-22 (JVP)
# ------------
#
# + TLMIN6, TLMAX6, TLMIN7, TLMAX7 re-calculated from $attGTIFilteredEvents and re-written in eventFile ($imgEvents)
#   Find attitude GTI file before use it and syntax errors corrected
#
# Version 2.57 - 2016-03-21 (PR, JVP)
# ------------
#
# + TLMIN6, TLMAX6, TLMIN7, TLMAX7 re-calculated from $attGTIFilteredEvents and re-written in eventFile ($finEvents)
#
# Version 2.56 - 2015-11-16 (PR, JVP)
# ------------
#
# + Flare background estimation in DET coordinates
#
# Version 2.55 - 2014-11-04 (JVP)
# ------------
#
# + Changes for SAS_14.0.0 integration (aligning pnEvents and 'epproc' SAS task)
#   - epreject : only for TI or BU (both withxrlcorrection = Y), but not for IM
#   - epevents : withrdpha = Y for IM and TI, but withrdpha = N for BU
#   - epfast   : back again but only for BU
#   - epnoise  : removed

# Version 2.54 - 2014-10-07 (JVP)
# ------------
#
# + forward PHA_RDCO after evlistcomb, change epframes parameters.
#   (aligning pnEvents and 'epchain' SAS task)
#
# Version 2.53 - 2014-10-03 (JVP)
# ------------
#
# + Changes for SAS_14.0.0 integration (aligning pnEvents and 'epchain' SAS task)
#   - epreject : withxrlcorrection = N for IMAGE, TIMING and BURST
#
# Version 2.52 - 2014-09-17 (JVP)
# ------------
#
# + Add 3 new parameters to 'epevents':
#       ctilongtermsoft = Y, withbackgroundgain = Y and withpatternoffset = Y
#       (informed by M. Smith through A. Ibarra)
#
# Version 2.51 - 2014-09-11 (JVP)
# ------------
#
# + Changes for SAS_14.0.0 integration:
#   - epreject : withxrlcorrection = Y for IMAGE, TIMING and BURST
#   - epevents : withrdpha = Y for IMAGE, TIMING and BURST
#   - epfast : removed.
#
# Version 2.50 - 2014-04-09 (LT)
# ------------
#
# + Fixed epreject, epevents, epfast to correctly apply energy scale corrections
#
# Version 2.49 - 2014-02-20 (EC)
# ------------
#
# + forward REJ_CORR after evlistcomb, change epframes parameters, PPS SCR #7176
#
# Version 2.48 - 2014-02-16 (EC)
# ------------
#
# + no epreject for burst mode, PPS SCR #7169
#
# Version 2.47 - 2014-01-22 (EC)
# ------------
#
# + loading corrections (modify epreject, add epfast), PPS SCR #7169
#
# Version 2.46 - 2013-08-26 (EC)
# ------------
#
# + RATE column from evselect must be non-zero for final epiclccorr call
#
# Version 2.45 - 2013-06-11 (EC)
# ------------
#
# + Merge changes from SOC version.
#
# version 2.44 - 2013-01-03 DLG
#
# - shorten intermediate SFBLC region file subset names to match MOS convention
#
# version 2.43 - 2012-12-03 DLG
#
# - merge faint and bright region files for generation of source-free background light curve
#
# version 2.422 - 2013-01-10 (EO)
#
# + Changed the parameter evenset in call to attcalc from $tmpEvents to $finEevents (line 957).
#
# version 2.421 - 2012-11-22 (EO)
# 
# - Added execution of command attcalc after execution of evselect (line 951)
#   if condition ($datamode eq "IMAGE.EL"),
#   to compute the full image of a mosaic observation.
# 
# version 2.42 - 2012-11-06 DLG
#
# - change references to $pn_timebinsize to $timebinsize
#
# version 2.41 - 2012-10-25 DLG
#
# - require valid nonzero flare GTI for final epiclcorr call
#
# version 2.40 - 2012-10-25 DLG
#
# - Modify plotting of SN output graph
#
# version 2.39 - 2012-10-24 DLG
#
# - evselect makeratecolumn=Y, disable tabcalc cal
#
# version 2.38 - 2012-10-23 DLG
#
# - added allcamera=Y to epiclccorr calls
#
# version 2.37 - 2012-10-17 DLG
#
# - check for existence of SFBLC region files in evselect
#
# version 2.36 - 2012-10-09 DLG
#
# - set writedss=Y on remaining SSP evselect calls
# = remove writedss param from epiclccorr
#
# version 2.35 - 2012-10-08 DLG
#
# - try writedss=Y on epiclccorr call
# + move tabcalc T_ELAPSED after epiclccorr
#
# version 2.34 - 2012-10-03 DLG
#
# - reenable epiclccorr, specify outset as temp file
#
# version 2.33 - 2012-10-01 DLG
#
# - disable epiclccorr calls pending GTI crash fix
#
# version 2.32 - 2012-09-27 DLG
#
# + correct epiclccorr param to 'applyabsolutecorrections'
#
# version 2.31 - 2012-09-26 DLG
#
# + Apply relative corrections to PN timeseries via epiclccorr()
#
# version 2.30 - 2012-09-26 DLG
#
# + undo 2.29 (moved to EPICSourceProducts)
# + correct evselect logic for region files
# + modify RATE column to real64
#
# version 2.29 - 2012-09-04 DLG
#
# + copy region file extensions to FBKTSR product
#
# version 2.28 - 2012-09-04 DLG
#
# + copy region file extensions to FBKTSR product
#
# version 2.27 - 2012-09-04 DLG
#
# + Add source regions to evselect for Optimised flare bkg time series
#
# version 2.26 - 2012-08-15 DLG
#
# + append SN_TO_BKGCUT table on background time series product
# 
# version 2.25 - 2012-07-27 DLG
#
# + modify SFBLC region code: source exclusion radius to be dependent on source counts
# 
# version 2.24 - 2012-07-23 DLG
#
# + disable 100cts limit on $rcut_pn
# 
# version 2.23 - 2012-07-18 DLG
#
# + add call to epnoise to flag noisy CCDs
# 
# version 2.22 - 2012-07-09 DLG
#
# + add ignorelowcnttail=yes parameter for bkgoptrate calls
# 
# version 2.21 - 2012-07-02 DLG
#
# + Refactor code for sparse image handling
#
# version 2.20 - 2012-06-25 DLG
#
# + Remove global attitude GTI filter constraint from in-band SF light curve
#
# version 2.19 - 2012-06-25 DLG
#
# + Add FLCUTTHR header keyword to opt flare light curve 
# 
# version 2.18 - 2012-06-22 DLG
#
#  + Do not perform source-free background analysis for sparse images < 1000 counts
# 
# version 2.17 - 2012-06-21 DLG
#
#  + Generate in-band source-free light curves FBKTSR as products (FITS, PDF)
# 
# version 2.16 - 2012-06-20 DLG
#
#  + Optimal cut threshold limited to 100 for flare GTI generation
# 
# version 2.16 - 2012-05-25 DLG
# 
#  + Add SN diagnostic for source-free bkgoptrate calc
# 
# version 2.15 - 2012-05-25 DLG
# 
#  + Add source-free flare background extraction
# 
# version 2.14 - 2012-04-26 DLG
# 
# + Merge epexposure call from saseval8.4
# 
# version 2.13 - 2011-03-07 DLG
# 
# + Add test call to bkgoptrate, start optimised flare background subtraction development
# 
# version 2.12 - 2009-10-05 DLG
#
# + Plot EPIC background flare light curves in log space (Mantis 52)
#
# version 2.11 - 2009-09-20 DLG
#
# + Added counting cycle report check, set $mingtisize accordingly
#
#
# version 2.10 - 2009-07-28 DLG
#
# + set mingtisize=0 for testing
#
# version 2.09a - 2009-06-19 DLG
# 
# + ignored if SAS_ATTITUDE set to "RAF" for slew processing
#
# version 2.09 - 2009-06-11 DLG
# ------------
#
# + epreject calls on imaging mode only per Mantis #46
#
#
# version 2.08 - 2009-05-15 DLG
# ------------
#
# + added support for epreject per Mantis feature req #35
#
#
# version 2.07 - 2009-01-12 DLG
# ------------
# 
# + Changed "withtempcorrection" to "Y" following advice from Hermann Brunner
#
#
# version 2.06 - 2006-07-10 DJF
# ------------
#
# + Fix for darck quadrant problem applied to OAL.  No longer need OBT_WARN test.
#   Replaced test with an echo of the OBT_WARN value to the log file.
#
# version 2.05 - 2006-04-04 DJF
# ------------
#
# + Test for OBT_WARN keyword in epframes output. Return ignore if it is found to be true.
#
# version 2.04 - 2005-11-16 IMS
# ------------
#
# + Incorporated DJF's shortened intermediate file names.
#
# version 2.03 - 2005-10-14 IMS
# ------------
#
# + TSTART obtained from flare time series header and subtracted from the times of the flare PDF plot.
# + Flare plot Y axis now RATE not COUNTS. Calculation of RATE is thus moved forward a bit, and all flare productions now moved within the test for existence of RATE extension.
#
# version 2.02 - 2005-08-12 RGW
# ------------
# + changes to flare plotting to work around 80-character filename limit
#   in 'fplot'
#
# version 2.01 - 2005-08-02 IMS
# ------------
#
# + PDF product now made from flare bkg time series.
# + Slightly changed event selection expression (events with PI < 150 eV are no longer kept; also the flag filter mask is given explicitly, not as XMMEA_EP).
#
# version 2.00 - 2005-05-09 DJF
# ------------
#
# + Add explicit perl module headers.  Previously these were supplied
#   at compile time.  This will make debugging and extending the modules
#   through additional perl libraries easier.
#
# + Code layout made more uniform with perltidy
# + Standerdized date format (ccyy-mm-dd)
# + Standerdized author list to use comma separator
# + Make use of perl $VERSION magic.  Now $Version = version = ''
##
# Version 1.69 - 16-Mar-2004 (DJF)
# ------------
#
# + Added dummy badpixset=temporary file to first call to badpixfind
#   to stop it creating it's default.
#
# Version 1.68 - 14-Jan-2004 (DJF)
# ------------
#
# + Use tabgtigen parameter mingtisize rather than evselect to exclude
#   small GTI from flare background
#
# Version 1.67 - 10-12-2003 (DJF)
# ------------
#
# + Adapted doCommand to use anonymous lists for list parameters
# + Fixed filter expression.  PI<-150 should have been PI<=150
#
# Version 1.66 - 05-Sept-2003 (DJF)
# -------------------------------
#  + Fixed backgroundrate in second call to badpixfind.  Was 0.0011 should have been 0.00011
#
# Version 1.65 - 11-Nov-2002 (DJF)
# -------------------------------
#
# + Updated evlistcomb columns to include OFFSETX
# + Added DLIMAP to evlistcomb othertables option
#
# Version 1.63 - 11-Nov-2002 (DJF)
# -------------------------------
#
# + Updated backgroundrate (0.0001 -> 0.0011) and hienergythresh (10.0 -> 12.0)
#   in second call to badpixfind
#
# Version 1.62 - 17-Oct-2002 (DJF)
# -------------------------------
#
# + Moved OFFSETS from othertables to maintable option
#
# Version 1.61 - 19-Feb-2002 (DH)
# -------------------------------
#
# + Add OFFSETS to the list of tables for evlistcomb to propagate.
#
# Version 1.60 - 15-Feb-2002 (DH)
# -------------------------------
#
# + Bug fixes for 1.60 .
#
# Version 1.59 - 7-Feb-2002 (DH)
# -------------------------------
#
# + Remove RAWY>10 cut on event lists in imaging mode.
#
# Version 1.58 - 6-Nov-2001 (DH)
# -------------------------------
#
# + New parameters for badpixfind for making of the background
#   mask. Values as recommended by Andy Read, via email from Michael
#   Freyburg dated 10/10/01.
#
# Version 1.57 - 9-Oct-2001 (DH)
# -------------------------------
#
# + Change energy thresholds in badpixfind for making of the background
#   mask.  Was 0.16/10.0, now 7.0/15.0 .
#
# Version 1.56 - 4-Oct-2001 (DH)
# -------------------------------
#
# + Fix bug that made flare screening only work for full frame modes.
# + Change threshold on number for good background pixels from 0 to 100.
#
# Version 1.55 - 2-Oct-2001 (DH)
# -------------------------------
#
# + Put back flare gti creation for small window mode.  Instead
#   inhibit flare gti creation if the number of suitable background
#   pixels found in a field is zero.
# + Increase flare gti cut from 6 to 10.
#
# Version 1.54 - 1-Oct-2001 (DH)
# -------------------------------
#
# + Scale integration time for flare histograms ($timebinsize)
#   to 2 times nominal for large window mode, and 35 times for
#   small window mode.
# + Do not create flare gti file for small window mode, though
#   a flare histogram is still created.
#
# Version 1.53 - 1-Oct-2001 (DH)
# -------------------------------
#
# + Replace dslatts call with the hasFITSExtension function.
#
# Version 1.52 - 2001-08-09 (DJF)
# ------------
#
# + Changed Flare light curve from intermediate to final product.
#
# Version 1.51 - 2001-06-15
# ------------
#
# + Change withtempcorrection parameter back to 'N'.
#
# Version 1.50 - 2001-06-06
# ------------
#
# + Change withtempcorrection=N parameter of epevents to 'Y'.
#
# Version 1.49 - 2001-06-05
# ------------
#
# + Re-work merging of HK and CCD specific gti data.  gtialign
#   was not being called properly.
#
# Version 1.48 - 2001-05-22
# ------------
#
# + Only apply RAWY>10 only for imaging mode data.
# + Add withtempcorrection=N parameter to epevents, in case the
#   default changes.
#
# Version 1.47 - 2001-05-16
# ------------
#
# + Bug fixes and refinements for changes made in version 1.46
#
# Version 1.46 - 2001-05-16 (DH)
# ------------
#
# + Remove unwanted columns from the EXPOSUnn extensions
#   of the pn event lists.
#
# Version 1.45 - 2001-04-17 (DH)
# ------------
#
# + Changes to main badpixfind call:
#   - Change loenergythresh from 0.16 to 0.14
#   - Change hicolthresh from 0.0015 to 0.00105
#
# Version 1.44 - 2001-03-16 (DH)
# ------------
#
# + Print out version number in performAction() for
#   tracking purposes.
#
# Version 1.43 - 2001-02-01 (DH)
# ------------
#
# + Change threshold for flare rejection rate cutoff from 36.4 to 6
#
# Version 1.42 - 2001-01-26 (DH)
# ------------
#
# + Restrict event energy range for flare rejection to 7-15 keV
#
# Version 1.41 - 2001-01-19 (DH)
# ------------
#
# + Change narrowerthanpsf parameter in standard call
#   to badpixfind from 1.25 to 1.5.
#
# Version 1.40 - 2001-01-19 (DH)
# ------------
#
# + hicolthesh values for the two calls to badpixfind
#   were not set correctly.  In particular, value for the
#   standard call was off by a factor of 10.  Values now
#   set as per Andy Read's email dated 13-Dec-2000.
#
# Version 1.39 - 2001-01-11 (DH)
# ------------
#
# + Remove GlobalHK from ignore rule.
#
# + Add InstrumentHK to start rule.
#
# Version 1.38 - 2000-12-19 (DH)
# ------------
#
# + Change hicolthesh in badpixfind masking back
#   to its nominal value.
#
# Version 1.37 - 2000-12-13 (DH)
# ------------
#
# + First production version.
#
# Declare list of instruments this module is interested in
@instList = qw(epn);

# Number of streams
sub numberOfStreams
{
    return numberOfExposures();
}

# Rules method
sub evaluateRules
{

    #
    return ignore()
      if ignored(
        module => 'ExpChooser'
        , instrument => thisInstrument
        , stream => thisStream
      );

    # ignore for slew processing
    return ignore() if ($ENV{'PCMS_ISSLEW'});

    #
    start()
      if complete(
        module => 'ExpChooser'
        , instrument => thisInstrument
        , stream => thisStream
      )
      and complete(
        module => 'InstrumentHK'
        , instrument => thisInstrument
        , stream => 1
      )
      and complete(
        module => 'GlobalHK'
        , instrument => 'all'
        , stream => 1
      );
}

# Action method
sub performAction
{
    info("Module version number: $version");
    my @extraopts;
    my @extraootopts;

    # Work out which exposure this stream maps on to
    my $exp_id = exposureID(
        instrument => thisInstrument
        , stream => thisStream
    );
    info("Exposure ID: $exp_id");

    #
    my $filter = getExposureProperty(
        exp_id => $exp_id
        , name => 'filter'
    );
    my $mode = getExposureProperty(
        exp_id => $exp_id
        , name => 'mode'
    );
    info("Filter: $filter");
    info("Instrument mode: $mode");

    # Find event files for this exposure
    my @evFiles = findFile(
        class => 'ODF'
        , instrument => thisInstrument
        , exp_id => $exp_id
        , content => '*events'
    );

    # Return immediately if there are no event lists
    return success() if $#evFiles < 0;

    # Find attitude HK file
    my $attFile = findFile(
        class => 'product'
        , instrument => 'all'
        , content => 'Attitude time series'
        , required => 'true'
    );

    # Find HK GTI file
    my $hkGTIFile = findFile(
        class => 'intermediate'
        , instrument => thisInstrument
        , content => 'HK GTI'
    );

    # Set mingtisize to 100 by default
    my $mingtisize = 100;

    # Check for counting cycle report
    # reset mingtisize if required
    my @ccxfiles = findFile(class => 'ODF'
			    ,instrument => thisInstrument
			    ,exp_id => $exp_id
			    ,content => 'Counting cycle report'
			    );
    
    unless (@ccxfiles) {
	info("No counting cycle report found. Setting mingtisize=0");
	$mingtisize = 0;
    }

    #
    my ( @finalEvents, @finalootEvents, @gtiFilter, @maskFilter, $datamode );
    my $goodPixels = 0;

    # Time binning to use when making background light curves
    #
    my $timebinsize = 10.0;

    # Loop through each event file
    foreach my $evFile0 (@evFiles)
    {

        # Find CCD and node numbers
        my $ccdid = readFITSKeyword(
            file => $evFile0
            , extension => 2
            , keyword => "CCDID"
        );
        my $quadrant = readFITSKeyword(
            file => $evFile0
            , extension => 2
            , keyword => "QUADRANT"
        );
        my $ccd = $quadrant * 3 + $ccdid + 1;
        $datamode = readFITSKeyword(
            file => $evFile0
            , extension => 2
            , keyword => "DATATYPE"
        );
#        my $ccdmode = readFITSKeyword(
#            file => $evFile0
#            , extension => 2
#            , keyword => "CCDMODE"
#        );
#        info("Processing events for CCD $ccd, mode $datamode.");
# Unused. IMS 2005-08-03
        my $evFile1 = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , ccd => $ccd
            , content => 'Processed events'
            , level => 1
        );
        my $gtiFile1 = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , ccd => $ccd
            , content => 'CCD GTI'
            , level => 1
        );

        # Run epframes
        #
        doCommand(
            'epframes'
            , set => $evFile0
            , eventset => $evFile1
            , gtiset => $gtiFile1
            , withsrcrawy => 'no'
          )
          or return exception();


		my %obt_warn = ( file =>  $evFile1
			, extension => 'EXPOSURE'
			, keyword => 'OBT_WARN'
		);
		if (hasFITSKeyword( %obt_warn ))
		{
			my $opt_warn = readFITSKeyword( %obt_warn );
			if ( $opt_warn =~ /$RE{boolean}{true}/ )
			{	
				info("OBT_WARN True");
			}
			else
			{	
				info("OBT_WARN False");
			}
		}
        # Run badpixfind to create a mask to be used for
        # generating a background lightcurve
        #
        my $dummyBadPixels = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , ccd => $ccd
            , content => 'Ignore This Bad pixels'
        );
        my $badPixMask = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , ccd => $ccd
            , content => 'Bright pixels mask'
        );
        doCommand(
            'badpixfind', eventset => $evFile1
            , badpixmap => $badPixMask
            , badpixset => $dummyBadPixels
            , withfovmask => 'Y'
            , threshabovebackground => 'Y'
            , hithresh => 0.00025
            , hicolthresh => 0.00007
            , narrowerthanpsf => 0
            , withbadpixmap => 'Y'
            , withsearchbadcolumn => 'Y'
            , loenergythresh => 7.0
            , hienergythresh => 15.0
            , backgroundrate => -1
            , thresholdlabel => 'rate'
          )
          or return exception();
        $goodPixels += readFITSKeyword(
            file => $badPixMask
            , extension => 1
            , keyword => "GOODPIX"
        );
        push( @maskFilter
            , "(CCDNR==${ccd} && MASK($badPixMask,0,0,RAWX,RAWY))" );

        # Run badpixfind
        #
        @extraopts = ();
        if ( $datamode eq "IMAGE.EL" )
        {
            my $newBadPixels = newFile(
                class => 'intermediate'
                , instrument => thisInstrument
                , exp_id => $exp_id
                , ccd => $ccd
                , content => 'Bad pixels'
            );
            doCommand(
                'badpixfind'
                , eventset => $evFile1
                , badpixset => $newBadPixels
                , hithresh => 0.0045
                , flickertimesteps => 1
                , thresholdlabel => 'rate'
                , withsearchbadcolumn => 'Y'
                , hicolthresh => 0.00105
                , backgroundrate => 0.00011
                , loenergythresh => 0.14
                , hienergythresh => 12.0
                , narrowerthanpsf => 1.5
              )
              or return exception();
            @extraopts = (
                getnewbadpix => 'Y'
                , badpixset => $newBadPixels
            );
        }

        # Flag bad pixels 
        # 
        doCommand(
            'badpix'
            , eventset => $evFile1
            , getuplnkbadpix => 'Y'
            , getotherbadpix => 'Y'
            , windowfilter => 'N'
            , @extraopts
          )
          or return exception();

	
#############
	# Run epreject only for TI or BU, but not for IMAGE - reduce low-energy PN noise
	if ( $datamode eq "IMAGE.EL" )
	{
	info("DEBUG: epreject task does not run for PN Imaging mode.");
	}
	elsif ( $datamode eq "TIME.EL" || $datamode eq "BURST.EL" )
	{
		# withxrlcorrection=Y for TIMING and BURST
		my $withxrlcorrection = 'Y';
		doCommand(
                      'epreject'
                      , eventset => $evFile1
                      , set => $evFile0
                      , sigma => 4.0
                      , withoffsetmap => 'Y'
		      , withxrlcorrection => $withxrlcorrection
                      )
		or return exception();
		
	info("DEBUG: Datamode $datamode found. epreject task appplied.");
	}
	
#	# epreject (withxrlcorrection => 'N') for IMAGE, TIMING and BURST
#            doCommand(
#                      'epreject'
#                      , eventset => $evFile1
#                      , set => $evFile0
#                      , sigma => 4.0
#                      , withoffsetmap => 'Y'
#		      , withxrlcorrection => 'N'
#                      )
#                or return exception();
#	    
#            # also run epnoise to reject soft X-ray noise in EPIC-pn
#	    
#	    # epnoise makes sense when it is passed 2 times.
#	    # Here it is run only 1 time => Check if it is needed or not in the future.
#	    # 'epnoise' removed on Oct/2014
#	    #
#	    if ( $datamode eq "IMAGE.EL" ) {
#	    doCommand(
#		      'epnoise'
#		      , set => $evFile1
#		      , identifynoisyframes => 'Y'
#		      , noisecut => 2
#		      , sigmacut => 3.0
#		      , savemasks => 'N'
#		      )
#		or return exception();
#	    }
#
        # Run epevents
        #
        my $evFile2 = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , ccd => $ccd
            , content => 'Processed events'
            , level => 2
        );

        # withrdpha=Y only for TIMING
        # my $withrdpha = 'N';
	#
        # if ($datamode eq "TIME.EL"){
	#
        #    $withrdpha = 'Y';
	#
	# }
	
	# withrdpha=Y for IMAGE, TIMING and BURST
	# my $withrdpha = 'Y';
	
        # withrdpha=N only for BURST
        my $withrdpha = 'Y';
	
         if ($datamode eq "BURST.EL"){
            $withrdpha = 'N';
	 }
#############	
	
	info("Starting epevents....");	
        doCommand(
            'epevents'
            , eventset => $evFile1
            , gainctiaccuracy => 2
            , reemissionthresh => 0
            , randomizeposition => 'Y'
            , randomizeenergy => 'Y'
            , withtempcorrection => 'Y'
            , withrdpha => $withrdpha
            , outset => $evFile2
	    , ctilongtermsoft => 'Y'
	    , withbackgroundgain => 'Y'
	    , withpatternoffset => 'Y'
	    , backgroundtres => 'Y'
	    , backgroundtbin => 100
          )
          or return exception();
	info("Completed epevents....");
        # Apply attitude correction
        #
        doCommand(
            'attcalc'
            , eventset => $evFile2
            , refpointlabel => 'pnt'
            , attitudelabel => 'ahf'
            , withmedianpnt => 'Y'
            , atthkset => $attFile
            , -w => 10
          )
          or return exception();

        # Filter events before final merging
        my $evFile3 = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , ccd => $ccd
            , content => 'Processed events'
            , level => 3
        );
#        my $expression = '#XMMEA_EP && (PI>150 || PI<=150)';
        my $expression = '(FLAG & 0x0fa0000)==0 && PI>150';
        doCommand(
            'evselect'
            , table => "${evFile2}:EVENTS"
            , withfilteredset => 'Y'
            , filteredset => $evFile3
            , destruct => 'Y'
            , expression => $expression
            , writedss => 'Y'
            , updateexposure => 'Y'
            , keepfilteroutput => 'Y'
          )
          or return exception();

        # Append to list of processed event files
        #
        push( @finalEvents, $evFile3 );

        # Accumulate expression for GTI filter
        #
        push( @gtiFilter, "(CCDNR==${ccd} && GTI($gtiFile1,TIME))" );

	if ($datamode eq "IMAGE.EL"){
	    info("Starting out-of-time epevents ....");
            my $evootFile2 = newFile(
                                     class => 'intermediate'
                                     , instrument => thisInstrument
                                     , exp_id => $exp_id
                                     , ccd => $ccd
                                     , content => 'Processed oot events'
                                     , level => 2
                                     );

            doCommand(
                'epevents'
                , eventset => $evFile1
                , gainctiaccuracy => 2
                , reemissionthresh => 0
                , randomizeposition => 'Y'
                , randomizeenergy => 'Y'
                , withtempcorrection => 'Y'
                , withrdpha => $withrdpha
                , outset => $evootFile2
                , ctilongtermsoft => 'Y'
                , withbackgroundgain => 'Y'
                , withpatternoffset => 'Y'
                , withoutoftime => 'Y'
                , backgroundtres => 'Y'
                , backgroundtbin => 100
              )
              or return exception();
	    info("Completed out-of-time epevents....");

            # Apply attitude correction
            #
            doCommand(
                'attcalc'
                , eventset => $evootFile2
                , refpointlabel => 'pnt'
                , attitudelabel => 'ahf'
                , withmedianpnt => 'Y'
                , atthkset => $attFile
                , -w => 10
              )
              or return exception();

            # Filter events before final merging
            my $evootFile3 = newFile(
                class => 'intermediate'
                , instrument => thisInstrument
                , exp_id => $exp_id
                , ccd => $ccd
                , content => 'Processed oot events'
                , level => 3
            );
#           my $expression = '#XMMEA_EP && (PI>150 || PI<=150)';
            my $expression = '(FLAG & 0x0fa0000)==0 && PI>150';
            doCommand(
                'evselect'
                , table => "${evootFile2}:EVENTS"
                , withfilteredset => 'Y'
                , filteredset => $evootFile3
                , destruct => 'Y'
                , expression => $expression
                , writedss => 'Y'
                , updateexposure => 'Y'
                , keepfilteroutput => 'Y'
              )
              or return exception();

            # Append to list of processed event files
            #
            push( @finalootEvents, $evootFile3 );

	}
    }

    # Combine per-CCD event lists into single listinfo("DEBUG: epreject task does not run for PN Imaging mode.");
    #
    @extraopts = ();
    my $tmpEvents;
    my $tmpootEvents;
    if ( $datamode eq "IMAGE.EL" )
    {
        $tmpEvents = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , content => 'Merged imaging events'
        );
	
       $tmpootEvents = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , content => 'Merged imaging oot events'
        );

        @extraopts = (
            imagingset => $tmpEvents
            , epnimgcolnames => [
                qw(TIME RAWX RAWY DETX DETY X Y PHA PI FLAG PATTERN PAT_ID PAT_SEQ OFFSETX)
              ]
            , epnimgcoltypes => [
                qw(double int16 int16 int16 int16 int32 int32 int16 int16 int32 int8 int16 int8 int16)
            ]
        );
        @extraootopts = (
            imagingset => $tmpootEvents
            , epnimgcolnames => [
                qw(TIME RAWX RAWY DETX DETY X Y PHA PI FLAG PATTERN PAT_ID PAT_SEQ OFFSETX)
              ]
            , epnimgcoltypes => [
                qw(double int16 int16 int16 int16 int32 int32 int16 int16 int32 int8 int16 int8 int16)
            ]
        );
	
	# evlistcomb and epexposure for OoT events only for Imaging (not for Timing or Burst)
	doCommand(
	    'evlistcomb'
	    , eventsets => [@finalootEvents]
	    , instrument => 'epn'
	    , maintable => [qw(EVENTS OFFSETS)]
	    , othertables => [qw(BADPIX EXPOSURE DLIMAP HKAUX)]
	    , mainattributes => 'REJ_CORR PHA_RDCO'
	    , @extraootopts
	    )
	or return exception();

	# + + + Add epexposure call w/time randomisation, as requested by Carlos (also for OoT events)
	doCommand(
	    'epexposure'
	    , eventsets => $tmpootEvents
	    , screenexposure => 'Y'
	    , spatialexposure => 'N'
	    , randomizetime => 'Y'
            )
        or return exception();

    }
    elsif ( $datamode eq "TIME.EL" || $datamode eq "BURST.EL" )
    {
        $tmpEvents = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , content => 'Merged timing events'
        );
        @extraopts = (
            timingset => $tmpEvents
            , epntimcolnames => [
                qw(TIME RAWX RAWY DETX DETY X Y PHA PI FLAG PATTERN PAT_ID PAT_SEQ OFFSETX)
              ]
            , epntimcoltypes => [
                qw(double int16 int16 int16 int16 int32 int32 int16 int16 int32 int8 int16 int8 int16)
            ]
        );
    }

    doCommand(
        'evlistcomb'
        , eventsets => [@finalEvents]
        , instrument => 'epn'
        , maintable => [qw(EVENTS OFFSETS)]
        , othertables => [qw(BADPIX EXPOSURE DLIMAP HKAUX)]
        , mainattributes => 'REJ_CORR PHA_RDCO'
        , @extraopts
      )
      or return exception();

    # + + + Add epexposure call w/time randomisation, as requested by Carlos

    doCommand(
	      'epexposure'
	      , eventsets => $tmpEvents
	      , screenexposure => 'Y'
	      , spatialexposure => 'N'
	      , randomizetime => 'Y'
	      ) or return exception();


    # Remove unwanted columns from EXPOSUnn extensions
    my @badColumns = (qw()); # disabled, as epexposure does this already?
    #my @badColumns = (qw(FRAMELIM NDSCLIN1 NDSCLIN2 NDSCLIN3 NDSCLIN4));
    my @badObjects;
    for ( my $ic = 1 ; $ic <= 12 ; $ic++ )
    {
        my $ext = sprintf 'EXPOSU%02d', $ic;
        info("Checking for extension $ext");
        if ( hasFITSExtension( file => $tmpEvents, extension => $ext ) )
        {
            my @columns = @badColumns;
            push @badObjects
              , map { $_ = "$tmpEvents:$ext:$_"; } @columns;
        }
    }
    if (@badObjects)
    {
        doCommand( 'dsrm', objects => [@badObjects] )
          or return exception();
    }

    # Fold in the HK gti files
    #
    my $hkGtiFileAligned = newFile(
        class => 'intermediate'
        , instrument => thisInstrument
        , exp_id => $exp_id
        , content => 'Aligned HK GTIs'
    );
    if ( fileExists( file => $hkGTIFile ) )
    {
        doCommand(
            "gtialign"
            , gtitable => $hkGTIFile . ":STDGTI"
            , eventset => $tmpEvents
            , outset => $hkGtiFileAligned
          )
          or return exception();
        foreach my $gtiFile (
            findFile(
                class => 'intermediate'
                , instrument => thisInstrument
                , exp_id => $exp_id
                , content => 'CCD GTI'
                , level => 1
            )
          )
        {
            my $fileInf
              = fileInfo( class => 'intermediate', name => $gtiFile );
            my $gtiExt = sprintf 'STDGTI%02d', $$fileInf{'ccd'};
            doCommand(
                "gtimerge"
                , tables => [ "$gtiFile", "$hkGtiFileAligned:$gtiExt" ]
                , mergemode => 'and'
                , withgtitable => 'N'
              )
              or return exception();
        }
    }

    # Filter event list on GTI
    #
    my $gtiExpr = join( " || ", @gtiFilter );
    my $finEvents;
    my $finootEvents;
    if ( $datamode eq "IMAGE.EL" )
    {
        $finEvents = newFile(
            class => 'product'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , content => 'EPIC PN imaging mode event list'
        );
        $finootEvents = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , content => 'EPIC PN imaging mode OOT event list'
        );
        doCommand(
                  'evselect'
                  , table => $tmpootEvents
                  , filteredset => $finootEvents
                  , withfilteredset => 'Y'
                  , expression => $gtiExpr
                  , keepfilteroutput => 'Y'
                  , writedss => 'Y'
                  , updateexposure => 'Y'
                  , destruct => 'true'
                  )
        or return exception();
    }
    elsif ( $datamode eq "TIME.EL" || $datamode eq "BURST.EL" )
    {
        $finEvents = newFile(
            class => 'product'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , content => 'EPIC timing mode event list'
        );
    }
    doCommand(
        'evselect'
        , table => $tmpEvents
        , filteredset => $finEvents
        , withfilteredset => 'Y'
        , expression => $gtiExpr
        , keepfilteroutput => 'Y'
        , writedss => 'Y'
        , updateexposure => 'Y'
        , destruct => 'true'
      )
      or return exception();
    
    # Added to compute full image from a mosaic observation
    
    if ($datamode eq "IMAGE.EL"){
	doCommand(
		'attcalc'
		,eventset      => $finEvents
		,refpointlabel => 'pnt'
		,attitudelabel => 'ahf'
		,withmedianpnt => 'Y'
		,atthkset      => $attFile
		,calctlmax     => 'yes'
		,-w            => 10
		)
		or return exception();

	doCommand(
		'attcalc'
		,eventset      => $finootEvents
		,refpointlabel => 'pnt'
		,attitudelabel => 'ahf'
		,withmedianpnt => 'Y'
		,atthkset      => $attFile
		,calctlmax     => 'yes'
		,-w            => 10
		)
		or return exception();

    ########################################################################
	my $attGTIFilteredEvents = newFile(
				class => 'intermediate'
				, instrument => thisInstrument
				, exp_id => $exp_id
				, content => 'Att GTI filtered events'
				);
	
	my $attGTIFilteredootEvents = newFile(
				class => 'intermediate'
				, instrument => thisInstrument
				, exp_id => $exp_id
				, content => 'Att GTI filtered oot events'
				);

	# Find attitude GTI file
	    my $attGTIFile = findFile(
				      class => 'intermediate'
				      , instrument => 'all'
				      , content => 'attitude GTI'
				      );
	doCommand(
		'evselect'
		, table => $finEvents
		, filteredset => $attGTIFilteredEvents
		, withfilteredset => 'Y'
		, expression => "GTI($attGTIFile,TIME)"
		)
		or return exception();

	doCommand(
		'attcalc'
		,eventset      => $attGTIFilteredEvents
		,refpointlabel => 'pnt'
		,attitudelabel => 'ahf'
		,withmedianpnt => 'Y'
		,atthkset      => $attFile
		,calctlmax     => 'yes'
		,-w            => 10
		)
		or return exception();

	doCommand(
		'evselect'
		, table => $finootEvents
		, filteredset => $attGTIFilteredootEvents
		, withfilteredset => 'Y'
		, expression => "GTI($attGTIFile,TIME)"
		)
		or return exception();

	doCommand(
		'attcalc'
		,eventset      => $attGTIFilteredootEvents
		,refpointlabel => 'pnt'
		,attitudelabel => 'ahf'
		,withmedianpnt => 'Y'
		,atthkset      => $attFile
		,calctlmax     => 'yes'
		,-w            => 10
		)
		or return exception();

	# Re-write TLMIN6, TLMAX6, TLMIN7, TLMAX7 keywords in the eventFile: $finEvents
	my $tlmin6_inter = readFITSKeyword(file => $attGTIFilteredEvents, extension => 'EVENTS', keyword => 'TLMIN6');
	my $tlmax6_inter = readFITSKeyword(file => $attGTIFilteredEvents, extension => 'EVENTS', keyword => 'TLMAX6');
	my $tlmin7_inter = readFITSKeyword(file => $attGTIFilteredEvents, extension => 'EVENTS', keyword => 'TLMIN7');
	my $tlmax7_inter = readFITSKeyword(file => $attGTIFilteredEvents, extension => 'EVENTS', keyword => 'TLMAX7');

	FITSUtils::writeKeyword($finEvents, "EVENTS", "TLMIN6" => {value => $tlmin6_inter});
	FITSUtils::writeKeyword($finEvents, "EVENTS", "TLMAX6" => {value => $tlmax6_inter});
	FITSUtils::writeKeyword($finEvents, "EVENTS", "TLMIN7" => {value => $tlmin7_inter});
	FITSUtils::writeKeyword($finEvents, "EVENTS", "TLMAX7" => {value => $tlmax7_inter});

	# Re-write TLMIN6, TLMAX6, TLMIN7, TLMAX7 keywords in the eventFile: $finootEvents
	FITSUtils::writeKeyword($finootEvents, "EVENTS", "TLMIN6" => {value => $tlmin6_inter});
	FITSUtils::writeKeyword($finootEvents, "EVENTS", "TLMAX6" => {value => $tlmax6_inter});
	FITSUtils::writeKeyword($finootEvents, "EVENTS", "TLMIN7" => {value => $tlmin7_inter});
	FITSUtils::writeKeyword($finootEvents, "EVENTS", "TLMAX7" => {value => $tlmax7_inter});
	#
	# We write down the same Legal Limits ( tlmin6/7 , tlmax6/7 ) in both $finEvents and $finootEvents
	# to ensure the same exactly borders for both event files.
	#
	# In principle the LegalLimits should be the same in both $attGTIFilteredEvents and $attGTIFilteredootEvents
	# but randomization could produce some mismatches.
	#
	# Be aware that TLMIN6, TLMAX6, TLMIN7, TLMAX7 are again re-written into the final products,
	# FWC and OOT imaging event files by the module 'EPICImageSize.pm'
	#

    ########################################################################


        # get FWC
        #
	# PN case : Search for FWC data only if PN in 'PrimeFullWindow' or 'PrimeFullWindowExtended' mode
	my $submode = readFITSKeyword(
				      file => $finEvents
				      , extension => 'PRIMARY'
				      , keyword => 'SUBMODE'
				      ) or return exception();

	if ( $submode eq "PrimeFullWindow" || $submode eq "PrimeFullWindowExtended" ) {

	        info("Imaging mode is PrimeFullWindow or PrimeFullWindowExtended. Producing FWC data...");

                my $finFWCEvents = newFile(
                    class => 'intermediate'
                    , instrument => thisInstrument
                    , exp_id => $exp_id
                    , content => 'EPIC PN imaging mode FWC event list'
                );
		doCommand(
		    'evqpb'
		    , table => $finEvents
		    , attfile => $attFile
		    , outset => $finFWCEvents
		    , exposurefactor => 4
		    , calctlmax => 'Y'
		    , overwritesubmode => 'Y'		# Full frame mode checking disabled. \
                                                        # This task should only be used in full frame images
		    )
		    or return exception();

		# Re-write TLMIN6, TLMAX6, TLMIN7, TLMAX7 keywords also to the FWC event list: $finFWCEvents
		FITSUtils::writeKeyword($finFWCEvents, "EVENTS", "TLMIN6" => {value => $tlmin6_inter});
		FITSUtils::writeKeyword($finFWCEvents, "EVENTS", "TLMAX6" => {value => $tlmax6_inter});
		FITSUtils::writeKeyword($finFWCEvents, "EVENTS", "TLMIN7" => {value => $tlmin7_inter});
		FITSUtils::writeKeyword($finFWCEvents, "EVENTS", "TLMAX7" => {value => $tlmax7_inter});

	} else {
	        info("$submode imaging mode is not Full Frame.");
		info("No FWC data produced for ", thisInstrument , ", exposure: $exp_id");
	}


    }





# epfast for burst modes
    if ( $datamode eq "BURST.EL" ) {
	doCommand(
		  'epfast'
		  , evfile => $finEvents
		  )
	    or return exception();
    }



    # Copy CIF into event list
    #
    my $cifFile = findFile(
        class => 'product'
        , content => 'Calibration index file'
    );
    doCommand(
        'dscopyblock'
        , blocks => "${cifFile}:CALINDEX"
        , to => $finEvents
      )
      or return exception();

    # Make flare GTI
    #
    if ( $datamode eq "IMAGE.EL" )
    {
        $timebinsize *= 2  if $mode eq "PrimeLargeWindow";
        $timebinsize *= 35 if $mode eq "PrimeSmallWindow";
        my $flareExpr = "($gtiExpr)";
        if (@maskFilter)
        {
            $flareExpr .= " && (" . join( " || ", @maskFilter ) . ")";
        }

        # Save copy of bright pixel/CCD GTI filter for SFBLC processing

	my $maskExpr = $flareExpr;	

        # Put in energy range for flare rejection

        $flareExpr .= " && (PI IN [7000:15000])";
        my $backRates = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , content => 'Flare background timeseries'
        );
        my $flareGti = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , content => 'flare GTI'
        );
        doCommand(
		  'evselect'
		  , table => $finEvents
		  , maketimecolumn => 'Y'
		  , makeratecolumn => 'Y'
		  , timebinsize => $timebinsize
		  , withrateset => 'Y'
		  , rateset => $backRates
		  , writedss => 'Y'
		  , updateexposure => 'N'
		  , expression => $flareExpr
		  )
	    or return exception();
	
#
# 1 pixel ~ 17 arcsec**2 ~ 0.00474 arcmin**2, convert to per [ks * arcmin**2]
#
        my $scalefactor = 0.0047355e-3;  # unit = [arcmin**2/pixel ks/s]

# If rate file is empty, no rate extension is created.  Trap this.
        info("Number of good background pixels found: $goodPixels");
        if ( hasFITSExtension( file => $backRates, extension => "RATE" )
            && $goodPixels > 100 )
        {
        #    doCommand(
        #        'tabcalc'
        #        , tables => "$backRates:RATE"
        #        , columntype => 'real64'
        #        , column => 'RATE'
        #        , expression =>
        #          "(COUNTS/($timebinsize*$goodPixels*$scalefactor))"
        #      )
        #      or return exception();

	    # Now subtract TSTART from the TIME values and leave the result in a new column T_ELAPSED:
            my $start_time=readFITSKeyword(
                file => $backRates
                , extension => 'RATE'
                , keyword => "TSTART"
              );
            doCommand(
                'tabcalc'
                , tables => "$backRates:RATE"
                , columntype => 'real64'
                , column => 'T_ELAPSED'
                , expression => "TIME-$start_time" # Use the perl variable rather than the string '#TSTART' in case we decide later to use some other time as offset.
              )
              or return exception();
            my $fplotCommandFile = newFile(class => 'intermediate'
	        , content => 'Flare bkg ts pco file'
                , format => 'ASCII'
	        , extension => 'pco'		      
            );
            writeASCIIFile(
                name => $fplotCommandFile
                , text => [join("\n", "LA X TIME-$start_time (s)", 'LA G2 RATE (cts/s)', 'log y', 'PLOT', 'QUIT', '')]
            );
            my $flareBkgPs = newFile(
                class => 'intermediate'
                , instrument => thisInstrument
                , exp_id => $exp_id
                , content => 'Flare bkg ts'
                , format => 'PS'
            );
            doCommand(
                'fplot'
                , infile => $backRates
                , xparm => 'T_ELAPSED'
                , yparm => 'RATE'
                , rows => '-'
                , device => "${flareBkgPs}/cps"
                , pltcmd => '@'."$fplotCommandFile"
              )
              or return exception();

#
# 1 pixel ~ 17 arcsec**2 ~ 0.00474 arcmin**2, convert to per [ks * arcmin**2]
#
#        my $scalefactor = 0.0047355e-3;  # unit = [arcmin**2/pixel ks/s]
#
# If rate file is empty, no rate extension is created.  Trap this.
#        info("Number for good background pixels found: $goodPixels");
#        if ( hasFITSExtension( file => $backRates, extension => "RATE" )
#            && $goodPixels > 100 )
#        {
#            doCommand(
#                'tabcalc'
#                , tables => "$backRates:RATE"
#                , columntype => 'real32'
#                , column => 'RATE'
#                , expression =>
#                  "(COUNTS/($timebinsize*$goodPixels*$scalefactor))"
#              )
#              or return exception();


# apply relative corrections to flare background timeseries intermediate

	    my $tmpBackRates = newFile(
				       class => 'intermediate'
				       , instrument => thisInstrument
				       , exp_id => $exp_id
				       , content => 'Copy flare background timeseries'
				       );
	    
	    doCommand(
		      'epiclccorr'
		      , srctslist => $backRates
		      , eventlist => $finEvents
		      , outset => $tmpBackRates
		      , applyabsolutecorrections => 'N'
		      , allcamera => 'Y'
		      ) or return exception();

	    copyFile(
		     source => $tmpBackRates
		     , destination => $backRates
		     );

# Add bkgoptrate test call, progress towards flare optimised background subtraction (Mantis 73)

	    my $retval = 0;
	    my $optCmd= "bkgoptrate fracexpstyle=auto tssettabname=$backRates:RATE withendtime=N withlowercutoffcount=N withmintimeratio=N ignorelowcnttail=Y withstarttime=N -w 10";
	    info("RUN: $optCmd");
	    my $optRate = `$optCmd`;
	    chomp $optRate;
	    $retval = $?;
	    if ($retval ne 0) { return exception("Error $retval returned by bkgoptrate: $!"); }
	    info("OPTRATE: $optRate");

            doCommand('tabgtigen', table => $backRates
		      , gtiset => $flareGti
		      , expression => "RATE<$optRate"
		      , mingtisize => $mingtisize
		      ) or return exception();

    # Insert source-free bkg flare light curve determination

# in-band source-free lightcurve generation for PN
# to include in pnEvents/sasflare8.5

	    info("==== Start in-band source-free flare lightcurve generation");

# Generate image for SFBLC analysis

	    my $SFBLCInputImage = newFile(
					  class => 'intermediate'
					  , instrument => thisInstrument
					  , exp_id => $exp_id
					  , content => 'SFBLC image'
					  , band => 8
					  , format => 'FITS'
					  );

# Find attitude GTI file
	    my $attGTIFile = findFile(
				      class => 'intermediate'
				      , instrument => 'all'
				      , content => 'attitude GTI'
				      );
	    
# define inband image filter expression
	    
	    my $inband_image_expr = "(RAWY > 12) && GTI($flareGti, TIME) && GTI($attGTIFile, TIME) && (( PATTERN <= 4 && PI in (500:1000] && (FLAG & 0x2fa002c) == 0) || (PATTERN<=4 && PI in (1000:7500] && (FLAG & 0x2fa0024) == 0))";
	    

# make image (make use of existing flare GTI)

	    # Define pix size
	    my $xbin_SFBLC = 80;
	    my $ybin_SFBLC = 80;

	    doCommand(
		      'evselect'
		      , table => $finEvents
		      , imageset => $SFBLCInputImage
		      , xcolumn => 'DETX'
		      , ycolumn => 'DETY'
		      , imagebinning => 'binSize'
		      , ximagebinsize => $xbin_SFBLC
		      , yimagebinsize => $ybin_SFBLC
		      , expression => $inband_image_expr
		      , imagedatatype => 'Int32'
		      , withimagedatatype => 'true'
		      , writedss => 'true'
		      , updateexposure => 'true'
		      , keepfilteroutput => 'false'
		      ) or return exception();


	    
#	    my $pn_timebinsize; downstream references to $pn_timebinsize changed to $timebinsize
	    my $sourceFreeBackRates;
	    my $sourceFreeImage;
	    my $SFBLCRegionFileFaint;
	    my $SFBLCRegionFileBright;
	    my $SFBLCBackScal;


	my $SFBLCDetmask;
	my $SFBLCDetMaskArith;
	my $source_free_expr;
	# Skip if input image $SFBLCInputImage has EXPOSURE = 0
	if ( readFITSKeyword(file => $SFBLCInputImage, extension => 'PRIMARY', keyword => 'EXPOSURE') != 0 ) {

#	    # check at least 1000 counts in input image before continuing to process
#           my $source_free_expr; # selection expression for source-free background light curve
#	    if (! &ModuleUtil::isImageLimit($SFBLCInputImage,'',1000.0)) {
#
#           The Mask is produced although image counts < 1000

# make vignetted exposure map
		
		my $SFBLCExpmap = newFile(
					  class => 'intermediate'
					  , instrument => thisInstrument
					  , exp_id => $exp_id
					  , content => 'SFBLC exposure map'
					  , band => 8
					  , format => 'FITS'
					  );
		
		doCommand(
			  'eexpmap'
			  , imageset => $SFBLCInputImage
			  , eventset => $finEvents
			  , attitudeset => $attFile
			  , attrebin => 2
			  , expimageset => $SFBLCExpmap
			  , pimin => 500
			  , pimax => 7500
			  , withvignetting => 'yes'
                          , withdetcoords => 'yes'
			  ) or return exception();
		
# make non-vignetted exposure map
		
		my $SFBLCNonVigExpmap = newFile(
						 class => 'intermediate'
						, instrument => thisInstrument
						, exp_id => $exp_id
						, content => 'SFBLC non-vig exposure map'
						, band => 8
						, format => 'FITS'
						);
		
		doCommand(
			  'eexpmap'
			  , imageset => $SFBLCInputImage
			  , eventset => $finEvents
			  , attitudeset => $attFile
			  , attrebin => 2
			  , expimageset => $SFBLCNonVigExpmap
			  , pimin => 500
			  , pimax => 7500
			  , withvignetting => 'no'
                          , withdetcoords => 'yes'
			  ) or return exception();
		
		
# make detection mask
		
		$SFBLCDetmask = newFile(
					   class => 'intermediate'
					   , instrument => thisInstrument
					   , exp_id => $exp_id
					   , content => 'SFBLC detmask'
					   , band => 8
					   , format => 'FITS'
					   );
		
		doCommand(
			  'emask'
			  , detmaskset => $SFBLCDetmask
			  , expimageset => $SFBLCNonVigExpmap
			  , threshold1 => 0.5
			  , threshold2 => 1
			  );

##########################################################################################################################################
			# check at least 1000 counts in input image before continuing to process: Source detection, etc

			#my $SFBLCDetMaskArith;

			#my $source_free_expr;
			if (! &ModuleUtil::isImageLimit($SFBLCInputImage,'',1000.0)) {

# now run eboxdetect in local mode
# not trying to merge images
		
		info("Running eboxdetect in local mode on EPIC PN, exposure $exp_id...");
		
# declare box-local source list
		
		my $SFBLCBoxList = newFile(
					   class => 'intermediate'
					    , instrument => thisInstrument
					   , exp_id => $exp_id
					   , content => 'SFBLC box-local source list'
					   , band => 8
					    , format => 'FITS'
					   );

# declare box-local source list (faint subset)

		my $SFBLCBoxListFaint = newFile(
						class => 'intermediate'
						, instrument => thisInstrument
						, exp_id => $exp_id
						, content => 'SFBLC box-local source list faint subset'
						, band => 8
						, format => 'FITS'
						);

# declare box-local source list (bright subset)

		my $SFBLCBoxListBright = newFile(
						 class => 'intermediate'
						 , instrument => thisInstrument
						 , exp_id => $exp_id
						 , content => 'SFBLC box-local source list bright subset'
						 , band => 8
						 , format => 'FITS'
						 );
###################################################################################################################################################		
# declare box-local source list (subset to remove (drop) in case of Moving Object = Solar System Object == SSO)

		my $sso_object = $ENV{PCMS_SSO_OBJECT};
		
		my $SFBLCBoxListDrop;
		if ( $sso_object ){
		info("DEBUG: Pipeline processing for a Moving Object = SSO. Object reference system. PCMS_SSO_OBJECT = $sso_object ");
		
		$SFBLCBoxListDrop = newFile(
						 class => 'intermediate'
						 , instrument => thisInstrument
						 , exp_id => $exp_id
						 , content => 'SFBLC box-local source list drop subset'
						 , band => 8
						 , format => 'FITS'
						 );
		}
		
# generate box-local source list
	    
		doCommand(
			  'eboxdetect'
			  , imagesets => $SFBLCInputImage
			  , expimagesets => $SFBLCNonVigExpmap
			  , detmasksets => $SFBLCDetmask
			  , boxlistset => $SFBLCBoxList
			  , usemap => 'false'
			  , likemin => 5
			  , boxsize => 5
			  , withdetmask => 'true'
			  , nruns => 1
			  , ecf => "1.0"
			  , pimin => 500
			  , pimax => 7500
			  ) or return exception();
		
# Now need to isolate ID_INT=ID_BAND=0 lines to a new file
# and convert to a region file

		doCommand(
			  'fselect'
			  , infile => $SFBLCBoxList
			  , outfile => $SFBLCBoxListFaint
			  , expr => "(RATE <= 0.35)"
			  ) or return exception();

		doCommand(
			  'fselect'
			  , infile => $SFBLCBoxList
			  , outfile => $SFBLCBoxListBright
			  , expr => "(RATE > 0.35)"
			  ) or return exception();
###################################################################################################################################################		
		if ( $sso_object ){
		doCommand(
			  'fselect'
			  , infile => $SFBLCBoxList
			  , outfile => $SFBLCBoxListDrop
			  , expr => "(SCTS > 30)"
			  ) or return exception();
		}
		
		my $SFBLCTempFile = newFile(
					    class => 'intermediate'
					    , instrument => thisInstrument
					    , exp_id => $exp_id
					    , content => 'SFBLC temp file'
					    , band => 8
					    , format => 'FITS'
					    );
		
		$SFBLCRegionFileFaint = newFile(
						class => 'intermediate'
						, instrument => thisInstrument
						, exp_id => $exp_id
						, content => 'SFBLC reg faint'
						, band => 8
						, format => 'FITS'
						);
		
		$SFBLCRegionFileBright = newFile(
						 class => 'intermediate'
						 , instrument => thisInstrument
						 , exp_id => $exp_id
						 , content => 'SFBLC reg bright'
						 , band => 8
						 , format => 'FITS'
						 );
###################################################################################################################################################		
		my $SFBLCRegionFileDrop;
		if ( $sso_object ){
		$SFBLCRegionFileDrop = newFile(
						 class => 'intermediate'
						 , instrument => thisInstrument
						 , exp_id => $exp_id
						 , content => 'SFBLC reg drop'
						 , band => 8
						 , format => 'FITS'
						 );
		}
		
		my $lthresh = 50; # ML threshold for removing bright sources  
		my $rad = 60; # Fixed extraction rafius in arcsec (faint)
		my $rad2 = 100; # Fixed extraction rafius in arcsec (bright)
		
		my $region_select = "(ID_BAND == 1 && ID_INST == 1 && LIKE >= $lthresh && BOX_CTS > 0)";
		my $region_select_drop = "(ID_BAND == 1 && ID_INST == 1 && LIKE >= 30 && BOX_CTS > 0)";
		
		doCommand(
			  'region'
			  , eventset => "${finEvents}:EVENTS"
			  , tempset => $SFBLCTempFile
			  , srclisttab => $SFBLCBoxListFaint
			  , expression => $region_select
			  , operationstyle => 'global'
			  , fixedradius => $rad
			  , bkgregionset => "${SFBLCRegionFileFaint}:REGION"
                          , outunit => "detxy"
			  ) or return exception();

		doCommand(
			  'region'
			  , eventset => "${finEvents}:EVENTS"
			  , tempset => $SFBLCTempFile
			  , srclisttab => $SFBLCBoxListBright
			  , expression => $region_select
			  , operationstyle => 'global'
			  , fixedradius => $rad2
			  , bkgregionset => "${SFBLCRegionFileBright}:REGION"
                          , outunit => "detxy"
			  ) or return exception();
###################################################################################################################################################			  
		if ( $sso_object ){
		doCommand(
			  'region'
			  , eventset => "${finEvents}:EVENTS"
			  , tempset => $SFBLCTempFile
			  , srclisttab => $SFBLCBoxListDrop
			  , expression => $region_select_drop
			  , operationstyle => 'global'
                          , radiusstyle => 'enfrac'
			  , bkgregionset => "${SFBLCRegionFileDrop}:REGION"
                          , outunit => "detxy"
			  ) or return exception();
		}
		
# Generate temporary evaluation image
		
		my $SFBLCEvalImage = newFile(
					     class => 'intermediate'
					     , instrument => thisInstrument
					     , exp_id => $exp_id
					     , content => 'SFBLC temp evaluation image'
					      , band => 8
					     , format => 'FITS'
					     );
		
# Include regions in eval selection expression

		my $eval_image_expr = $inband_image_expr . " && (region($SFBLCRegionFileFaint,DETX,DETY) && region($SFBLCRegionFileBright,DETX,DETY))";
		
		doCommand(
			  'evselect'
			  , table => $finEvents
			  , imageset => $SFBLCEvalImage
			  , xcolumn => 'DETX'
			  , ycolumn => 'DETY'
			  , withimageset => 'Y'
			  , imagebinning => 'binSize'
			  , ximagebinsize => $xbin_SFBLC
			  , yimagebinsize => $ybin_SFBLC
			  , expression => $eval_image_expr
			  , imagedatatype => 'Int32'
			  , withimagedatatype => 'true'
			  , writedss => 'true'
			  , updateexposure => 'true'
			  , keepfilteroutput => 'false'
			  ) or return exception();
		
# Now, finally should be able to use the region file on the event list to 
# generate the light curve for flare filtering
		
# PN
# need to establish the exposure ID and other parts of the name
# needs the CCD GTI and bright pixel masks
# Am following the PN light curve creation steps which look more
# complicated/extensive than for MOS
	    
#	    info("Creating CCD/GTI filters for EPIC PN");
		
#	    my @gtiexpr;
		
#	    foreach my $gtiFile (
#				 findFile(
#					  class => 'intermediate'
#					  , instrument => thisInstrument
#					  , exp_id => $exp_id
#					  , content => 'CCD GTI'
#					  , level => 1
#					  )
#				 )
#	    {
#		my $fileInf
#		    = fileInfo( class => 'intermediate', name => $gtiFile );
#		my $ccd = $$fileInf{'ccd'};
#		push  @gtiexpr, "( CCDNR==$ccd && GTI($gtiFile, TIME))";
#	    };
	    
#	    my $gtiexpr_string = "(" . join(" || ", @gtiexpr) . ")";
	    
# 	    
# 	    info("Creating Bright Pixel Masks for EPIC PN");
# 	    
# 	    my @brtexpr;
# 	    
# 	    foreach my $brtmask (
# 				 findFile(
# 					  class => 'intermediate'
# 					  , instrument => thisInstrument
# 					  , exp_id => $exp_id
# 					  , content => 'Bright pixels mask'
# 					  )
# 				 )
# 	    {
# 
# 		my $fileInf
# 		    = fileInfo( class => 'intermediate', name => $brtmask );
# 		my $ccd = $$fileInf{'ccd'};
# 		push @brtexpr, "( CCDNR == $ccd && MASK($brtmask,1,1,RAWX,RAWY))";
# 	    };
	    
# 	    my $brtexpr_string = "(" . join(" || ", @brtexpr) . ")";

# merge faint and bright region files prior to generating source free light curve

		my $SFBLCRegionFileMerged = newFile(
						    class => 'intermediate'
						    , instrument => thisInstrument
						    , exp_id => $exp_id
						    , content => 'SFBLC reg merged'
						    , band => 8
						    , format => 'FITS'
						    );

		doCommand(
			  'fmerge'
			  , infiles => [$SFBLCRegionFileBright, $SFBLCRegionFileFaint]
			  , outfile => $SFBLCRegionFileMerged
			  , columns => "SHAPE,DETX,DETY,R,ROTANG,COMPONENT"
			  , mextname => 'REGION'
			  , history => 'Y'
			  , copyprime => 'Y'
			  ) or return exception();

		doCommand(
			  'fparkey'
			  , value => 'pos'
			  , fitsfile => "${SFBLCRegionFileMerged}[REGION]"
			  , keyword => "MTYPE1"
			  , add => 'Y'
			  ) or return exception();

		doCommand(
			  'fparkey'
			  , value => 'DETX,DETY'
			  , fitsfile => "${SFBLCRegionFileMerged}[REGION]"
			  , keyword => "MFORM1"
			  , add => 'Y'
			  ) or return exception();


# Calculate AREASCALE parameter in background estimate for flaring filtering
		my $SFBLCNonVigExpmapRegMask = newFile(
						 class => 'intermediate'
						, instrument => thisInstrument
						, exp_id => $exp_id
						, content => 'SFBLC non-vig exposure map reg mask'
						, band => 8
						, format => 'FITS'
						);
		$SFBLCDetMaskArith = newFile(
						 class => 'intermediate'
						, instrument => thisInstrument
						, exp_id => $exp_id
						, content => 'SFBLC det mask arith'
						, band => 8
						, format => 'FITS'
						);


		my $SFBLCRegionFileMerged_rows = readFITSKeyword(
								file => $SFBLCRegionFileMerged
								, extension => 'REGION'
								, keyword => 'NAXIS2'
								);
		if ( $SFBLCRegionFileMerged_rows ){
			info("DEBUG - $SFBLCRegionFileMerged has $SFBLCRegionFileMerged_rows rows");
		doCommand(
			  'regionmask'
			  , region => "${SFBLCRegionFileMerged}[REGION]"
			  , pixdefset => $SFBLCNonVigExpmap
			  , maskset => $SFBLCNonVigExpmapRegMask
			  , whichpixdef => 'image'
			  , withmaskset => 'Y'
			  ) or return exception();

		### farith ~/unam.fits intermediate/pnEvents-epn-1-SFBLC_detmask-S00300000080000.ftz\[MASK\] MUL ~/otram.fits
		doCommand('farith'
			  ,infil1 => $SFBLCNonVigExpmapRegMask
			  ,infil2 => "${SFBLCDetmask}[MASK]"
			  ,outfil => $SFBLCDetMaskArith
			  ,ops => 'MUL'
			  ) or return exception();
		} else {
			info("DEBUG - $SFBLCRegionFileMerged has 0 rows");
			$SFBLCDetMaskArith = "${SFBLCDetmask}[MASK]";
			#info("DEBUG - SFBLCDetMaskArith = $SFBLCDetMaskArith");
		}


# make in-band source-free light curve

		info("Creating source-free in band light curve for PN");

#		$source_free_expr = $maskExpr." && (PATTERN <=4) && (PI IN [500:7500]) && ((FLAG & 0xfffffef) == 0) && (region($SFBLCRegionFileFaint,X,Y) && region($SFBLCRegionFileBright,X,Y))";
		$source_free_expr = $maskExpr." && (PATTERN <=4) && (PI IN [500:7500]) && ((FLAG & 0xfffffef) == 0) && (region($SFBLCRegionFileMerged,DETX,DETY))";

	    } else {
		 # extraction expression for sparse images (omit region select)

		info("Creating in band light curve for PN (sparse image < 1000 counts)");

		$source_free_expr = $maskExpr." && (PATTERN <=4) && (PI IN [500:7500]) && ((FLAG & 0xfffffef) == 0)";

		#  No region file for sources removing in sparse image < 1000 counts
		info("DEBUG - No sources removing in sparse image < 1000 counts");
		$SFBLCDetMaskArith = "${SFBLCDetmask}[MASK]";

	    }
	} else {
	        info("Input SFBLC image $SFBLCInputImage has EXPOSURE = 0");
	        info("Creating in band light curve for PN (EXPOSURE = 0)");
		$source_free_expr = $maskExpr." && (PATTERN <=4) && (PI IN [500:7500]) && ((FLAG & 0xfffffef) == 0)";
	}


		if ( $SFBLCDetMaskArith ){
		#if (fileExists( file => $SFBLCDetMaskArith ) ){
		#
		#  Calculating SFBLCBackScal from $SFBLCDetMaskArith
		#	Cases:
		#	- Image counts > 1000
		#		- With Region file	: $SFBLCDetMaskArith == ${SFBLCDetmask}[MASK] * $SFBLCNonVigExpmapRegMask
		#		- With no Region file	: $SFBLCDetMaskArith == ${SFBLCDetmask}[MASK]
		#
		#	- Sparse image counts < 1000
		#		- With no Region file	: $SFBLCDetMaskArith == ${SFBLCDetmask}[MASK]
		#
		### fimgstat ~/otram.fits INDEF INDEF
		my %SFBLCDetMaskStat = &ModuleUtil::getImageStats($SFBLCDetMaskArith,'');
			my $xbin = readFITSKeyword(
						file => $SFBLCDetmask
						, extension => 'MASK'
						, keyword => 'CDELT1L'
						);
			my $ybin = readFITSKeyword(
						file => $SFBLCDetmask
						, extension => 'MASK'
						, keyword => 'CDELT2L'
						);

		$SFBLCBackScal = $SFBLCDetMaskStat{sum} * $xbin * $ybin;
		info("DEBUG - SFBLCBackScal = $SFBLCBackScal");
		#foreach my $stat_param ( sort { $a <=> $b } keys(%SFBLCDetMaskStat) ){
		#    info("DEBUG - SFBLCDetMaskStat $stat_param = $SFBLCDetMaskStat{$stat_param}");
		#}
		}


	    # Generate source-free background timeseries as product

#	    $pn_timebinsize = 10;
		
	    $sourceFreeBackRates  = newFile(
					    class => 'product'
					    , instrument => thisInstrument
					    , exp_id => $exp_id
					    , content => 'EPIC flare background timeseries' 
					    );
 	    $sourceFreeImage  = newFile(
					    class => 'intermediate'
					    , instrument => thisInstrument
					    , exp_id => $exp_id
					    , content => 'SFBLC source free image' 
					    );

	    doCommand(
		      'evselect'
		      , expression => $source_free_expr
		      , maketimecolumn => 'Y'
		      , makeratecolumn => 'Y'
		      , rateset => $sourceFreeBackRates
		      , table => $finEvents
		      , timebinsize => $timebinsize
		      , updateexposure => 'N'
		      , withrateset => 'Y'
		      , writedss => 'Y'
                      , withimageset => 'Y'
                      , ximagebinsize => $xbin_SFBLC
                      , yimagebinsize => $ybin_SFBLC
                      , xcolumn => 'DETX'
                      , ycolumn => 'DETY'
                      , imagebinning => 'binSize' 
                      , imageset => $sourceFreeImage
		      ) or return exception();

# add light curve rate and T_ELAPSED columns. The RATE column definition
# differs from current pipeline value because there the cut is made by
# rate per unit sky area. Here we don't care about unit area so we can
# just divide the counts by the bin time

# 	    info("Adding column for RATE, using PN time binsize $pn_timebinsize");
# 	    
# 	    doCommand(
# 		      'tabcalc'
# 		      , column => 'RATE'
# 		      , columntype => 'real64'
# 		      , expression => "COUNTS/${pn_timebinsize}"
# 		      , tables => "${sourceFreeBackRates}:RATE"
# 		      ) or return exception();
	    
# Observation start time should already be determined in pnEvents
	    
# Apply relative corrections to output product timeseries

	    my $tmpSourceFreeBackRates = newFile(
						 class => 'intermediate'
						 , instrument => thisInstrument
						 , exp_id => $exp_id
						 , content => 'Copy EPIC flare background timeseries' 
						 );

	    doCommand(
		      'epiclccorr'
		      , srctslist => $sourceFreeBackRates
		      , eventlist => $finEvents
		      , outset => $tmpSourceFreeBackRates
		      , applyabsolutecorrections => 'N'
		      , allcamera => 'Y'
		      ) or return exception();

	    copyFile(
		     source => $tmpSourceFreeBackRates
		     , destination => $sourceFreeBackRates
		     );


# Add BACKSCAL (AREASCALE) keyword to RATE header extension
		    if ( defined($SFBLCBackScal) ){
		    doCommand(
			      'addattribute'
			      , set => "${sourceFreeBackRates}:RATE"
			      , attributename => 'BACKSCAL'
			      , attributetype => 'real'
			      , attributecomment => '"\"Backscale for flaring filtering\""'
			      , realvalue => $SFBLCBackScal
			      ) or return exception();
		    }

# Generate T_ELAPSED column

	    info("Adding column for T_ELAPSED, using observation start time $start_time");
	    
	    doCommand(
		      'tabcalc'
		      , column => 'T_ELAPSED'
		      , columntype => 'real64'
		      , expression => "TIME-$start_time"
		      , tables => "${sourceFreeBackRates}:RATE"
		      ) or return exception();

# Run bkgoptrate - optimum cut threshold

	    my $snDiagFile = newFile(
				     class => 'intermediate'
				     , instrument => thisInstrument
				     , exp_id => $exp_id
				     , content => 'Flare GTI SN diagnostic'
				     , band => 8
				     , format => 'FITS'
				     );

	    my $retval = 0;
	    my $rcut_pn = `bkgoptrate -V 0 -w 0 tssettabname=${sourceFreeBackRates}:RATE dooutputsntab=yes ignorelowcnttail=yes snsettabname=${snDiagFile}:SN_TO_BKGCUT`;
	    info("CMD: bkgoptrate -V 0 -w 0 tssettabname=${sourceFreeBackRates}:RATE 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_pn eq '') {
		info("bkgoptrate call returns no value, skipping optimised flare GTI generation...");
	    } else {

		if (fileExists( file => $sourceFreeBackRates ) ) {

		    info("bkgoptrate returns optimum cut threshold: $rcut_pn");
		    
		    # append SN_TO_BKGCUT table to FBKTSR product
		    
		    doCommand(
			      'ftappend'
			      , infile => "${snDiagFile}[SN_TO_BKGCUT]"
			      , outfile => $sourceFreeBackRates
			      ) or return exception();
		    
		    # add FLCUTTHR header value
		    
		    doCommand(
			      'addattribute'
			      , set => "${sourceFreeBackRates}:RATE"
			      , attributename => 'FLCUTTHR'
			      , attributetype => 'real'
			      , attributecomment => '"\"Optimised flare cut threshold\""'
			      , realvalue => $rcut_pn
			      ) or return exception();
		    
# cut threshold limited to 100
		    
#		if ($rcut_pn > 100) {
#		    $rcut_pn = 100;
#		    info("Optimum cut threshold limited to 100");
#		}
		    
# Now generate the optimised flare GTI files
		    
		    $optFlareGti = newFile(
					   class => 'intermediate'
					   , instrument => thisInstrument
					   , exp_id => $exp_id
					   , content => "Optimised flare GTI"
					   , band => 8
					   , format => 'FITS'
					   );
		    
		    info("Creating GTI files with this threshold");
		    
		    doCommand(
			      'tabgtigen'
			      , table => $sourceFreeBackRates
			      , gtiset => $optFlareGti
			      , expression => "RATE < ${rcut_pn}"
			      , mingtisize => 100
			      ) or return exception();
	    
		}

	    }

# make plot of background light curve
	    
	    info("Make full in-band, source-free light curve plots");
	    #
            if ( hasFITSExtension(file => $sourceFreeBackRates, extension => "RATE") &&
                 hasFITSExtension(file => $sourceFreeBackRates, extension => "SN_TO_BKGCUT") ){
                 my $flareBkgPDF0 = newFile(
                                      class => 'product'
                                      , instrument => thisInstrument
                                      , exp_id => $exp_id
                                      , content => 'EPIC flare background timeseries'
                                      , format => 'PDF'
                                      );
                 my $flareBkgPNG0 = newFile(
                                      class => 'product'
                                      , instrument => thisInstrument
                                      , exp_id => $exp_id
                                      , content => 'EPIC flare background timeseries'
                                      , format => 'PNG'
                                      );
                 doCommand(
                           'plotFlrBkgTimeSeries.py'
                           , rateset  => $sourceFreeBackRates
                           , plotfile => $flareBkgPNG0
                           , pdfplotfile => $flareBkgPDF0
                           , plotformat => "both"
                           ) or return exception();
            } else {
                    info("No RATE or SN_TO_BKGCUT extensions found in source-free light curve. No plot created.");
            }


# Regenerate light curve but filter by new GTIs and plot to see if they disappear

	    info("TESTING: Recreating light curves with only valid GTIs and regions left in");
	    
	    my $optBackRates = newFile(
				       class => 'intermediate'
				       , instrument => thisInstrument
				       , exp_id => $exp_id
				       , content => 'Optimised flare bkg time series'
				       , band => 8
				       , format => 'FITS'
				       );

	    my $test_ts_expr = $maskExpr." && (PATTERN <=4) && GTI($attGTIFile, TIME) && (PI in [500:7500]) && ((FLAG & 0xfffffef) == 0)";

	    if (fileExists( file => $SFBLCRegionFileFaint )) {
		$test_ts_expr .= " && region($SFBLCRegionFileFaint,DETX,DETY)";
	    }
	    
	    if (fileExists( file => $SFBLCRegionFileBright )) {
		$test_ts_expr .= " && region($SFBLCRegionFileBright,DETX,DETY)";
	    }
	    
	    if (fileExists( file => $optFlareGti )) {
		$test_ts_expr .=  " && GTI($optFlareGti, TIME)";
	    }

	    doCommand(
		      'evselect'
		      , expression => $test_ts_expr
		      , maketimecolumn => "Y"
		      , makeratecolumn => 'Y'
		      , rateset => $optBackRates
		      , table => $finEvents
		      , timebinsize => $timebinsize
		      , updateexposure => 'N'
		      , withrateset => 'Y'
		      , writedss => 'Y'
		      ) or return exception();
	    
# 	    doCommand(
# 		      'tabcalc'
# 		      , column => 'RATE'
# 		      , columntype => 'real64'
# 		      , expression => "COUNTS/${pn_timebinsize}"
# 		      , tables => "${optBackRates}:RATE"
# 		      ) or return exception();
	    
# Apply relative corrections to test GTI timeseries

	    my $tmpOptBackRates = newFile(
					  class => 'intermediate'
					  , instrument => thisInstrument
					  , exp_id => $exp_id
					  , content => 'Copy optimised flare bkg time series'
					  , band => 8
					  , format => 'FITS'
					  );

	    # Require valid nonzero flare GTI, and non-zero RATE column  for final epiclccorr call

	    if( (numberFITSRows(file=>$optFlareGti, extension=>'STDGTI') > 0) &&
                !ModuleUtil::isColumnEmpty($optBackRates,"RATE") )
	    {

		doCommand(
			  'epiclccorr'
			  , srctslist => $optBackRates
			  , eventlist => $finEvents
			  , outset => $tmpOptBackRates
			  , applyabsolutecorrections => 'N'
			  , allcamera => 'Y'
			  ) or return exception();
		
		
		copyFile(
			 source => $tmpOptBackRates
			 , destination => $optBackRates
			 );
		
		doCommand(
			  'tabcalc'
			  , column => 'T_ELAPSED'
			  , columntype => 'real64'
			  , expression => "TIME-$start_time"
			  , tables => "${optBackRates}:RATE"
			  ) or return exception();


                # add FLCUTTHR header value. Needed for plotFlrBkgTimeSeries.py
                doCommand(
                              'addattribute'
                              , set => "${optBackRates}:RATE"
                              , attributename => 'FLCUTTHR'
                              , attributetype => 'real'
                              , attributecomment => '"\"Optimised flare cut threshold\""'
                              , realvalue => $rcut_pn
                              ) or return exception();

                my $OptSourceFreeLightCurvePDF = newFile(
                                                        class => 'intermediate'
                                                        , instrument => thisInstrument
                                                        , exp_id => $exp_id
                                                        , content => 'Opt in-band source-free light curve'
                                                        , band => 8
                                                        , format => 'PDF'
                                                        );
		doCommand(
                           'plotFlrBkgTimeSeries.py'
                           , rateset  => $optBackRates
                           , pdfplotfile => $OptSourceFreeLightCurvePDF
                           , plotformat => "pdf"
                           ) or return exception();

	    }
        }
    }

    # Everything OK
    return success();
}
1;