Download

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.46';
$version = $VERSION;
$author  = 'Dean Hinshaw,Richard West,Duncan John Fyfe,Ian Stewart,Duncan Law-Green,Ed Chapin';
$date    = '2013-08-26';

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

    # 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, @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
            , wrongpixlimit => 20
            , srcposition => 190
            , eventset => $evFile1
            , gtiset => $gtiFile1
            , mipmethod => 'onboard'
            , mipthreshold => 3000
            , setupbpx => 'nom3'
          )
          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 - reduce low-energy PN noise
        # Imaging mode only

	if ( $datamode eq "IMAGE.EL" ) {
	    doCommand(
		      'epreject'
		      , eventset => $evFile1
		      , set => $evFile0
		      , sigma => 4.0
#		      , noiseparameter => "0.98 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0"
		      , withoffsetmap => 'Y'
		      )
		or return exception();

        # run epnoise to reject soft X-ray noise in EPIC-pn
	    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
        );
	info("Starting epevents....");
        doCommand(
            'epevents'
            , eventset => $evFile1
            , gainctiaccuracy => 2
            , reemissionthresh => 0
            , randomizeposition => 'Y'
            , randomizeenergy => 'Y'
            , withtempcorrection => 'Y'
            , outset => $evFile2
          )
          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))" );
    }

    # Combine per-CCD event lists into single list
    #
    @extraopts = ();
    my $tmpEvents;
    if ( $datamode eq "IMAGE.EL" )
    {
        $tmpEvents = newFile(
            class => 'intermediate'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , content => 'Merged imaging 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)
            ]
        );
    }
    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)]
        , @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;
    if ( $datamode eq "IMAGE.EL" )
    {
        $finEvents = newFile(
            class => 'product'
            , instrument => thisInstrument
            , exp_id => $exp_id
            , content => 'EPIC PN imaging mode event list'
        );
    }
    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();

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

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

# check at least 1000 counts in input image before continuing to process

	    my $source_free_expr; # selection expression for source-free background light curve 
	    
#	    my $pn_timebinsize; downstream references to $pn_timebinsize changed to $timebinsize
	    my $sourceFreeBackRates;
	    my $SFBLCRegionFileFaint;
	    my $SFBLCRegionFileBright;

	    if (! &ModuleUtil::isImageLimit($SFBLCInputImage,'',1000.0)) {

# 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'
			  ) 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'
			  ) or return exception();
		
		
# make detection mask
		
		my $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
			  );
		
# 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'
						 );
		
# 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();
		
		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 $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)";
		
		doCommand(
			  'region'
			  , eventset => "${finEvents}:EVENTS"
			  , tempset => $SFBLCTempFile
			  , srclisttab => $SFBLCBoxListFaint
			  , expression => $region_select
			  , operationstyle => 'global'
			  , fixedradius => $rad
			  , bkgregionset => "${SFBLCRegionFileFaint}:REGION"
			  ) or return exception();

		doCommand(
			  'region'
			  , eventset => "${finEvents}:EVENTS"
			  , tempset => $SFBLCTempFile
			  , srclisttab => $SFBLCBoxListBright
			  , expression => $region_select
			  , operationstyle => 'global'
			  , fixedradius => $rad2
			  , bkgregionset => "${SFBLCRegionFileBright}:REGION"
			  ) 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,X,Y) && region($SFBLCRegionFileBright,X,Y))";
		
		doCommand(
			  'evselect'
			  , table => $finEvents
			  , imageset => $SFBLCEvalImage
			  , xcolumn => 'X'
			  , ycolumn => 'Y'
			  , withimageset => 'Y'
			  , imagebinning => 'binSize'
			  , ximagebinsize => 80
			  , yimagebinsize => 80
			  , 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,X,Y,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 => 'X,Y'
			  , fitsfile => "${SFBLCRegionFileMerged}[REGION]"
			  , keyword => "MFORM1"
			  , add => 'Y'
			  ) or return exception();
	    
# 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,X,Y))";

	    } 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)";
		
	    }

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

	    doCommand(
		      'evselect'
		      , expression => $source_free_expr
		      , maketimecolumn => 'Y'
		      , makeratecolumn => 'Y'
		      , rateset => $sourceFreeBackRates
		      , table => $finEvents
		      , timebinsize => $timebinsize
		      , updateexposure => 'N'
		      , withrateset => 'Y'
		      , writedss => 'Y'
		      ) 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
		     );

# 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");
	    
	    my $sourceFreeLightCurvePS = newFile(
						 class => 'intermediate'
						 , instrument => thisInstrument
						 , exp_id => $exp_id
						 , content => "In-band source-free light curve"
						 , band => 8
						 , format => 'PS'
						 );
	    
	    my $sourceFreeLCCommandFile = newFile(
						  class => 'intermediate'
						  , instrument => thisInstrument
						  , exp_id => $exp_id
						  , content => 'In-band sf lc pco file'
						  , format => 'ASCII'
						  , extension => 'pco'
						  );

            # Cutoff threshold label string
	    my $srcut = sprintf("%.3f",$rcut_pn);
	    
            # Write PLT command file
	    writeASCIIFile(
			   name => $sourceFreeLCCommandFile
			   , text => [join("\n","LA X TIME-$start_time (s)",'LA G2 RATE (cts/s)','log y','lwid 5',"LA 1 P 0 $srcut LI 0 1 CO 4 \"$srcut\"",'PLOT','QUIT','')]
			   );
	    
	    doCommand(
		      'fplot'
		      , infile => $sourceFreeBackRates
		      , xparm => 'T_ELAPSED'
		      , yparm => 'RATE'
		      , rows => '-'
		      , device => "${sourceFreeLightCurvePS}/CPS"
		      , pltcmd => '@'."$sourceFreeLCCommandFile"
		      ) or return exception();

# Get extremal values of calculated S/N

	    my $snBg = readFITSColumn(
				      file => $snDiagFile
				      , extension => 'SN_TO_BKGCUT'
				      , column => 'SN_RATIO'
				      ) or return exception();
	    
	    my $snMin = min(@$snBg);
	    my $snMax = max(@$snBg);
	    
# Generate plot of S/N

	    my $SNRCurvePS = newFile(
				     class => 'intermediate'
				     , instrument => thisInstrument
				     , exp_id => $exp_id
				     , content => 'In-band sf Signal to noise plot'
				     , band => 8
				     , format => 'PS'
				     );
	    
	    my $SNRCommandFile = newFile(
					 class => 'intermediate'
					 , instrument => thisInstrument
					 , exp_id => $exp_id
					 , content => 'In-band sf Signal to noise pco file'
					 , format => 'ASCII'
					 , extension => 'pco'
					 );
	    
	    writeASCIIFile(
			   name => $SNRCommandFile
			   ,text => [join("\n","LA X Background rate cut (cts/s)",'LA G2 Signal-to-noise','lwid 5',"LA 1 Pos $srcut $snMax LI -90 1.0 CO 4 \"$srcut\"",'PLOT','QUIT','')]
			   );
	    
	    doCommand(
		      'fplot'
		      , infile => $snDiagFile
		      , xparm => 'BKGRATECUT'
		      , yparm => 'SN_RATIO'
		      , rows => '-'
		      , device => "${SNRCurvePS}/CPS"
		      , pltcmd => '@'."$SNRCommandFile"
		      ) or return exception();

	    my $tempPS = newFile(
				 class => 'intermediate'
				 , instrument => thisInstrument
				 , exp_id => $exp_id
				 , content => 'In-band sf temp plot file'
				 , band => 8
				 , format => 'PS'
				 );
	    
	    my $combiPS = newFile(
				  class => 'intermediate'
				  , instrument => thisInstrument
				  , exp_id => $exp_id
				  , content => 'In-band sf combi plot file'
				  , band => 8
				  , format => 'PS'
				  );
	    
	    $retval=0;
	    my $cli_1 = `gs -sDEVICE=pswrite -sOutputFile=${tempPS} -dNOPAUSE -dBATCH $sourceFreeLightCurvePS $SNRCurvePS`;
	    info("CMD: gs -sDEVICE=pswrite -sOutputFile=${tempPS} -dNOPAUSE -dBATCH $sourceFreeLightCurvePS $SNRCurvePS");
	    $retval = $?;
	    if ($retval ne 0) { return exception("Error returned by gs: $!"); }
	    
	    $retval=0;
	    my $cli_2 = `psnup -s0.6 -2up ${tempPS} > ${combiPS}`;
	    info("CMD: psnup -s0.6 -2up ${tempPS} > ${combiPS}");
	    $retval=$?;
	    if ($retval ne 0) { return exception("Error returned by psnup: $!"); }
	    
	    # Generate optimised flare background PDF
	    
            my $flareBkgPdf = newFile(
				      class => 'product'
				      , instrument => thisInstrument
				      , exp_id => $exp_id
				      , content => 'EPIC flare background timeseries'
				      , format => 'PDF'
				      );
            PStoPDF(
		    source => $combiPS
		    , destination => $flareBkgPdf
		    )
		or return exception();
	    
# 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,X,Y)";
	    }
	    
	    if (fileExists( file => $SFBLCRegionFileBright )) {
		$test_ts_expr .= " && region($SFBLCRegionFileBright,X,Y)";
	    }
	    
	    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();
		
		my $OptSourceFreeLightCurvePS = newFile(
							class => 'intermediate'
							, instrument => thisInstrument
							, exp_id => $exp_id
							, content => 'Opt in-band source-free light curve'
							, band => 8
							, format => 'PS'
							);
		
		doCommand(
			  'fplot'
			  , infile => $optBackRates
			  , xparm => 'T_ELAPSED'
			  , yparm => 'RATE'
			  , rows => '-'
			  , device => "${OptSourceFreeLightCurvePS}/CPS"
			  , pltcmd => '@'."$sourceFreeLCCommandFile"
			  ) or return exception();
		
	    }
        }
    }

    # Everything OK
    return success();
}
1;