package MOSEvents;
use strict;
use English;
use Carp;
use List::Util qw(min max);
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='3.55';
$version = $VERSION;
$author='Dean Hinshaw,Richard West,Duncan John Fyfe,Ian Stewart,Duncan Law-Green,Eduardo Ojero,Ed Chapin,Jose V Perea';
$date='2018-07-25';

#
# ChangeLog
# =========
#
# Version 3.55 - 2018-07-25 (JVP)
# ------------
#
# + Writing down the same Legal Limits ( tlmin6/7 , tlmax6/7 ) also in the FWC imaging event list: $imgFWCEvents
#   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 3.54 - 2018-06-14 (JVP)
# ------------
#
# + Search for FWC data (evqpb) always in MOS case because peripheral CCDs are always in 'PrimeFullWindow' mode.
#
# Version 3.53 - 2018-06-05 (JVP)
# ------------
#
# + Search for FWC data (evqpb) only if the observing mode is 'PrimeFullWindow'. No 'PrimeFullWindowExtended' mode for MOS
#
# Version 3.52 - 2017-11-24 (JVP)
# ------------
#
# + Check the input TS for 'bkgoptrate'. If all the RATE column values are 0 or constant, do not run 'bkgoptrate'.
#   Several parts of the code have to be skipped if the previoues condition occurs.
#   New subroutine created to check the RATE column values from a TS: 'check_min_max_rate_col'
#
# Version 3.51 - 2017-01-24 (JVP)
# ------------
#
# + The Mask is produced although image counts < 1000.
#   SFBLCBackScal calculated although image counts < 1000. Calculated only from the Mask.
#
# Version 3.50 - 2016-12-16 (JVP)
# ------------
#
# + Calculate AREASCALE parameter in background estimate for flaring filtering
#   New SAS task included: regionmask
#   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 3.49 - 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 3.48 - 2016-03-21 (PR, JVP)
# ------------
#
# + TLMIN6, TLMAX6, TLMIN7, TLMAX7 re-calculated from $attGTIFilteredEvents and re-written in eventFile ($imgEvents)
#
# Version 3.47 - 2015-11-16 (PR, JVP)
# ------------
#
# + Flare background estimation in DET coordinates
#
# Version 3.46 - 2014-06-25 (JVP)
# ------------
#
# + Plot for 'Flare GTI SN diagnostic' produced from extension [SN_TO_BKGCUT] of product 'EPIC flare background timeseries'
#   Filename labels in plots are final products now, instead intermediate files.
#
# Version 3.45 - 2013-08-26 (EC)
# ------------
#
# + RATE column from evselect must be non-zero for final epiclccorr call
#
# Version 3.44 - 2013-06-11 (EC)
# ------------
#
# + Merge changes from SOC version.
#
# version 3.431 - 2012-11-22 EO
# ----------------------------
#
# - Added execution of command attcalc after execution of evselect (line 1225)
#   if condition (fileExists(file=>$tmpimgEvents) is fulfilled (line 1217),
#   to compute the full image of a mosaic observation.
#
# version 3.43 - 2012-11-05 DLG
# -----------------------------
#
# - merge faint and bright region files for generation of source-free background light curve
#
# version 3.42 - 2012-11-05 DLG
# -----------------------------
#
# - require valid flare GTI for final epiclccorr call (crashes on empty GTI)
#
# version 3.41 - 2012-10-25 DLG
# -----------------------------
#
# - Modify plotting of SN output graph
#
# version 3.40 - 2012-10-24 DLG
# -----------------------------
#
# - evselect makeratecolumn=Y, disable subsequent tabcalc
#
# version 3.39 - 2012-10-22 DLG
# -----------------------------
#
# - all allcamera = Y to epiclccorr call
#
# version 3.38 - 2012-10-17 DLG
# -----------------------------
#
# - check for existence of SFBLC region files in evselect
#
# version 3.37 - 2012-10-09 DLG
# -----------------------------
#
# - set writedss=Y on remaining SSP evselect calls
# - remove writedss param from epiclccorr
#
# version 3.36 - 2012-10-08 DLG
# -----------------------------
#
# - try writedss=Y on epiclccorr call
# + move tablcalc T_ELAPSED after epiclccorr
#
# version 3.35 - 2012-10-01 DLG
# -----------------------------
#
# - reenable epiclccorr, specify outset as temp file
#
# version 3.34 - 2012-10-01 DLG
# -----------------------------
#
# - disable epiclccorr call pending GTI crash fix
#
# version 3.33 - 2012-09-27 DLG
# -----------------------------
#
# + correct epiclccorr param to 'applyabsolutecorrections'
#
# version 3.32 - 2012-09-26 DLG
# -----------------------------
#
# + Apply relative time series corrections using epiclccorr()
#
# version 3.31 - 2012-09-26 DLG
# -----------------------------
#
# + Undo change 3.30 (moved to EPICSourceProducts)
# + correct evselect logic for region files
# + modify RATE column to real64
#
# version 3.30 - 2012-09-04 DLG
# -----------------------------
#
# + copy region file extensions to FBKTSR product
#
# version 3.29 - 2012-09-04 DLG
# -----------------------------
#
# + Alter time binning on flare background timeseries to 26s (from 52s)
#
# version 3.28 - 2012-09-04 DLG
# -----------------------------
#
# + Add source regions to evselect for Optimised flare bkg time series
# 
# version 3.27 - 2012-08-22 DLG
# -----------------------------
#
# + Generate SFBLC and SNR dual plot, with optimised SN marked
# 
# version 3.26 - 2012-08-15 DLG
# -----------------------------
#
# + append SN_TO_BKGCUT table on background time series product
# 
# version 3.25 - 2012-07-27 DLG
# -----------------------------
#
# + modify SFBLC region code: source exclusion radius to be dependent on source counts
# 
# version 3.24 - 2012-07-23 DLG
# -----------------------------
#
# + disable 100cts maximum on $rcut_mos
# 
# version 3.23 - 2012-07-10 DLG
# -----------------------------
#
# + reinstate bkgoptrate call for bad pixel filtering
# 
# version 3.22 - 2012-07-09 DLG
# -----------------------------
#
# + add ignorelowcnttail=yes parameter for bkgoptrate calls
# 
# Version 3.21 - 2011-07-02 DLG
# -----------------------------
# 
# + Refactor code for sparse image handling
# 
# Version 3.20 - 2011-07-02 DLG
# -----------------------------
# 
# + Trap null return from bkgoptrate, skip optimised GTI generation
# 
# Version 3.19 - 2011-06-25 DLG
# -----------------------------
# 
# + Remove global attitude GTI filter constraint fron in-band SF light curve
# 
# Version 3.18 - 2011-06-25 DLG
# -----------------------------
# 
# + Add FLCUTTHR header keyword to opt flare light curve
# 
# Version 3.17 - 2011-06-22 DLG
# -----------------------------
# 
# + Do not perform source-free background analysis for sparse images < 1000 counts
# 
# Version 3.16 - 2011-06-21 DLG
# -----------------------------
# 
# + Produce optimised flare background time series FITS and PDF as product (FBKTSR)
# 
# Version 3.15 - 2011-06-20 DLG
# -----------------------------
# 
# + Clamp maximum threshold to 100 for flare GTI generation
# 
# Version 3.14 - 2011-06-10 DLG
# -----------------------------
# 
#  + Disable bkgoptrate call for bad pixel threshold
# 
# Version 3.13 - 2011-06-07 DLG
# -----------------------------
# 
#  + Add SN diagnostic for source-free bkgoptrate calc
#  
# Version 3.12 - 2011-05-25 DLG
# -----------------------------
# 
#  + Add source-free background flare lightcurve generation
# 
# Version 3.11 - 2011-04-06 DLG
# -----------------------------
# 
# + Add second bkgoptrate test, 1st for bright pixel check, 2nd for optimised GTI
# 
# Version 3.10 - 2011-03-07 DLG
# ------------
#  
# + Add bkgoptrate test, optimised flare background filtering (Mantis 73)
# 
# Version 3.091 - 2013-01-09 EC
# -------------
#
# + Add calls to manualintervention to allow skipping of CCDs
#
# Version 3.09 - 2009-10-05 DLG
# ------------
#
# + Plot EPIC background flare light curves in log space (Mantis 52)
#
# Version 3.08 - 2009-09-20 DLG
# ------------
# 
# + Added counting cycle report check, set $mingtisize accordingly
#
# Version 3.07 - 2009-07-28 DLG
# ------------
# 
# + Set mingtisize = 0 for testing
#
# Version 3.06a - 2009-06-19 DLG
# -------------
# 
# + Set ignore() if SAS_ATTITUDE environment variable is 'RAF' for slew processing
#
# Version 3.06 - 2005-11-16 IMS
# ------------
#
# + Incorporated DJF's shortened intermediate file names.
#
# Version 3.05 - 2005-10-31 IMS
# ------------
#
# + Parameter --getotherbadpix of badpix call changed from 'no' to 'yes' on advice from JB (should get rid of bright pixels in ccd 4, mos1.)
#
# Version 3.04 - 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.
# + Added 'PStoPDF', 'writeASCIIFile' to the list of methods called from ModuleResources.
#
# Version 3.03 - 2005-09-30 RGW
# ------------
#
# + OFFSETS table is propagated by evlistcomb
#
# version 3.02 - 2005-08-12 RGW
# ------------
#
# + trimmed filenames in flare plotting to work around FPLOT filename length
#   restriction
#
# version 3.01 - 2005-08-09 IMS
# ------------
#
# + PDF product now made from flare bkg time series.
# + Calls to fmedian removed because embadpixfind>2.0 now calculates this internally. Parameter --medianset of embadpixfind has likewise been removed.
#
# version 3.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.
#
# Version 2.33 - 2004-03-19 (DJF)
# ------------
#
# + Added 'FLICKERING' and 'ON_BADOFFSET' to flareFlags and eventlistFlags.
#
# Version 2.32 - 2004-03-17 (DJF)
# ------------
#
# + Added Explicit use of withimageset=>'N' for evselect.
#
# Version 2.31 - 2004-02-12 (DJF)
# ------------
#
# + Use new PCMS function gtiTimeThreshold() to return threshold
#
# Version 2.30 - 2004-01-14 (DJF)
# ------------
#
# + Use tabgitgen parameter mingtisize rather than evselect to remove 
#   short intervals from bad pixel and flare background GTI.
#
# Version 2.29 - 2003-12-09 (DJF)
# ------------
#
# + Adapted doCommand to use anonymous lists for list parameters
#
# Version 2.28 - 2003-07-01 (DJF)
# ------------
#
# + #XMMEA_EM selection had been applied to wrong evselect.  This caused
#   removal of GATTI events before  flare background construction.
#   Try again.  #XMMEA_EM selection applied to final event list
#
# Version 2.27 - 2003-07-01 (DJF)
# ------------
#
# + XMMEA_EM selection still not working, more testing.
#
# Version 2.26 - 2003-07-01 (DJF)
# ------------
#
# + Fixed use of #XMMEA_EM flag selection in good event filtering.
#   From (FLAG & #XMMEA_EM)==0 to (#XMMEA_EM)==0
#
# Version 2.25 - 2003-07-01 (DJF)
# ------------
#
# + Corrected good event filtering.  Had missed out the ==0 from the flag filter.
#
# Version 2.24 - 2003-06-30 (DJF)
# ------------
#
# + Added emosdatamodes => 'IMAGING' to first evlistcomb call.
#
# Version 2.23 - 2003-04-11 (DJF)
# ------------
#
# + Fixed command line to evlistcomb (instrument= 'emos' rather than thisInstrument)
#
# Version 2.22 - 2003-04-08 (DJF)
# ------------
#
# + Added Jean Ballet's newer recipy for improved flare screening
#
# Version 2.21 - 2002-07-01 (DJF)
# ------------
#
# + Changed image creation flags to XMMEA_EM
#
# Version 2.20 - 2002-03-21 (DH)
# ------------
#
# + Bug fix for version 2.19.
#
# Version 2.19 - 2002-03-20 (DH)
# ------------
#
# + Put in removal of short (<100s) flare GTIs for bad pixel
#   finding.
# + Use offset variance files in first call to emevents.
# + Add bitmask fitlering to final event list, so it makes it in
#   the data subspace.
#
# Version 2.18 - 2002-03-13 (DH)
# ------------
#
# + Fix bugs in identifying and handling of CCDs in low gain mode.
#
# Version 2.17 - 2002-03-11 (DH)
# ------------
#
# + Bug fix: was using wrong file name for output flare gti file.
#
# Version 2.16 - 2002-03-05 (DH)
# ------------
#
# + Fix error in background rate calculation.  Do rate calculation
#   based on the ccd area available, rather than the number of CCDs.
#
# Version 2.15 - 2002-02-20 (DH)
# ------------
#
# + Back out change in 2.14, as the MOS software is not yet ready to
#   handle this.
#
# Version 2.14 - 2002-02-19 (DH)
# ------------
#
# + Add OFFSETS to the list of tables for evlistcomb to propagate.
#
# Version 2.13 - 2002-02-18 (DH)
# ------------
#
# + Explicitly list all parameters for the fmedian ftool call.
# + More bug fixes
#
# Version 2.12 - 2002-02-15 (DH)
# ------------
#
# + Bug fixes for 2.11 .
#
# Version 2.11 - 2002-02-07 (DH)
# ------------
#
# + New method for specifying bitmask filters for event lists, in order
#   to make it more clear which flags are being used.  
# + Add in the UNDERSHOOT flag, for event list and flare screening.
# + The flare background timeseries is now generated as a rates file, not
#   a histogram.  The RATE column is normalized to 1/[ks*acrmin**2], and 
#   the GTI cut is made at 2.0 .
# + Check for GAIN_CCD keyword in frame file.  If set to 'LOW', do not
#   process events for this ccd.
#
# Version 2.10 - 2002-01-08 (DH)
# ------------
#
# + Test for XMMEA_22 in event list before making flare histograms.
#
# Version 2.09 - 2001-11-14 (DH)
# ------------
#
# + Fix bug in flare GTI application introduced in 2.07.
#
# Version 2.08 - 2001-11-06 (DH)
# ------------
#
# + Slight change to flare GTI test to more exactly mimic MakeImage.pm .
#
# Version 2.07 - 2001-10-31 (DH)
# ------------
#
# + Test flare GTI as in MakeImage.pm to see if it is above a time threshold before using it.
#
# Version 2.06 - 2001-08-20 (DJF)
# ------------
#
# + Added test for empty GTI when filtering out flares before searching for bright pixels.
#
# Version 2.05 - 2001-08-15 (DJF)
# ------------
#
# + Added test for column 'HISTO' to first call to tabgtigen. 
#
# Version 2.04 - 2001-08-14 (DJF)
# ------------
#
# + Remove flagtruncatede1=N flag from first call to emevents (to fix XMMEA_22 failure fo evselect)
# + Add flag findbright=N to first call to embadpixfind as suggesteg by Jean Ballet. 
#
# Version 2.03 - 2001-08-13 (DJF)
# ------------
#
# + Another fix for histogram intermadiate to product change.
#
# Version 2.02 - 2001-08-13 (DJF)
# ------------
#
# + Fix to flare histogram intermadiate to product change.
#
#
# Version 2.01 - 2001-08-09 (DJF)
# ------------
#
# + Brought up to date with changes to previous versions of MOSEvents.
# + Changed Flare Histogram from intermediate to final product.
# - Removed supurfluous findFile at end of module v1.13
#
#############################################
##  Changes Applied to v1.10 before     #####
##  implimentation of v 2.0.1           #####
##
## Version 1.13 - 2001-06-07 (DH)
## ------------
##
## + Take out changes in version 1.12 and move test for
##   empty event list to MakeImage module.
##
## Version 1.12 - 2001-06-07 (DH)
## ------------
##
## + Set ignore flag if final imaging event list is empty.
##
## Version 1.11 - 2001-06-07 (DH)
## ------------
##
## + Do not make flare gti if flare histogram is empty.
##
############################################
#
# Version 2.00 - 2001-05-29 (DH)
# ------------
#
# + Complete re-work of bad pixel finding using the new embadpixfind package.
#
# Version 1.10 - 2001-05-21 (DH)
# ------------
#
# + Changes made as instructed by Jean Ballet:
#   - changes to rejection of e3 and e4 events due to task changes
#   - Add call to emenergy before badpixfind in Imaging mode only
#   - Change to emenergy options for TIME and CTIME modes
#
# Version 1.09 - 2001-03-16 (DH)
# ------------
#
# + Print out version number in performAction() for
#   tracking purposes.
#
# Version 1.08 - 2001-01-11 (DH)
# ------------
#
# + Remove GlobalHK from ignore rule.
#
# + Add InstrumentHK to start rule.
#
# Version 1.07 - 2000-11-23 (DH)
# ------------
#
# + First production version.
#

# Declare list of instruments this module is interested in
@instList=qw(emos1 emos2);

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

sub evaluateRules {
    
    # If upstream module has been ignored, so if this one
    return ignore() if ignored(module => 'ExpChooser',
			       instrument => thisInstrument,
			       stream => thisStream);

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

    # Start conditions
    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);
}

my $flareFlags =   ['UNDERSHOOT'
	, 'BAD_E3E4'
	, 'ON_BADROW'
	, 'OUTSIDE_THRESHOLDS'
	, 'OUT_OF_CCD_WINDOW'
	, 'ON_BADPIX'
	, 'COSMIC_RAY'
	, 'IN_BAD_FRAME'
	, 'OUT_OF_FOV'
	, 'FLICKERING'
	, 'ON_BADOFFSET'
];

my $eventlistFlags =['UNDERSHOOT'
	, 'BAD_E3E4'
	, 'ON_BADROW'
	, 'OUTSIDE_THRESHOLDS'
	, 'OUT_OF_CCD_WINDOW'
	, 'ON_BADPIX'
	, 'COSMIC_RAY'
	, 'IN_BAD_FRAME'
	, 'FLICKERING'
	, 'ON_BADOFFSET'
];

sub performAction {
	info("Module version number: $version");
	# Misc. variables
	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 files
	if ($#evFiles<0)
	{
		ppd("000005");
		return success();
	}
	my @evFiles2;

	# Find auxiliary file for this exposure in the ODF
	my $auxFile=findFile(class => 'ODF'
		,instrument => thisInstrument
		,exp_id => $exp_id
		,content => 'Auxiliary data'
		,required => 'true'
	);

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

	# Loop through each event file
	my @filter;
	my $fovArea = 0;

    foreach my $evFile0 (@evFiles) 
	{

		# Find CCD and node numbers
		my $ccd=readFITSKeyword(file => $evFile0,
			extension => 2,
			keyword => "CCDID"
		);

                # manualintervention will check for keyword/keyval
                # pairs provided to the "--ignore" command-line option
                # of "sequence" that match the supplied "instrument"
                # and "exp_id" values, respectively. Although
                # manualintervention was originally provided to skip
                # entire exposures, there is no reason we can't use it
                # here, provided that the keyword for skipping CCDs of
                # an instrument is not identical to a pure instrument
                # name which would cause things to be skipped earlier
                # in ExpChooser. To skip CCDs we append "ccd" to the
                # instrument name.
                #   e.g.,
                #           sequence --ignore="emos1ccd=4,7"
                #   will skip CCDs 4 and 7 of MOS1.

                if (manualintervention( instrument => thisInstrument . "ccd", exposure => $ccd ))
                {
                    info("Manual Intervention. Ignoring CCD $ccd");
                    next;
                }


		my $node=readFITSKeyword(file => $evFile0,
			extension => 2,
			keyword => "CCDNODE"
		);
		my $datamode=readFITSKeyword(file => $evFile0,
			extension => 2,
			keyword => "DATATYPE"
		);
		info("First pass of events for CCD $ccd, node $node, mode=$datamode");

		# Run emframes
		#
		my $frameFile=newFile(class => 'intermediate',
			instrument => thisInstrument,
			exp_id => $exp_id,
			node => $node, 'ccd' => $ccd,
			content => 'Rectified frames'
	 	);
		my $gtiFile1=newFile(class => 'intermediate',
			instrument => thisInstrument,
			exp_id => $exp_id,
			node => $node,
			ccd => $ccd,
			content => 'CCD GTI',
			level => 1
		);
		my $evFile1=newFile(class => 'intermediate',
			instrument => thisInstrument,
			exp_id => $exp_id,
			node => $node,
			ccd => $ccd,
			content => 'Processed events',
			level => 1
		);
		@extraopts=();
		@extraopts=(flagbadtimes => 'yes', ingtiset => $hkGTIFile)
			if $hkGTIFile
		;

		doCommand('emframes',
			  auxiliaryset => $auxFile,
			  frameset => $frameFile,
			  withodfeventset => 'Y',
			  odfeventset => $evFile0,
			  neweventset => 'Y',
			  outeventset => $evFile1,
			  writegtiset => 'Y',
			  outgtiset => $gtiFile1,
			  @extraopts
		) or return exception();

		if( hasFITSKeyword(file=>$frameFile, extension=>'FRAMES', keyword=>'GAIN_CCD')
		   && readFITSKeyword(file=>$frameFile, extension=>'FRAMES', keyword=>'GAIN_CCD') =~ /^LOW/ 
		){
			info("CCD $ccd is in low gain mode.  Will not process events for this ccd.");
			ppd(id => "000007" , ccd => $ccd , node => $node );
			next;
		}

		push @evFiles2, $evFile0;

			#
		@extraopts=();
		if( $datamode eq "IMAGE.EL" || $datamode eq "RIMAGE.EL" )
		{
			my $evFile2=newFile(class => 'intermediate',
				instrument => thisInstrument,
				exp_id => $exp_id,
				node => $node,
				ccd => $ccd,
				content => 'Processed events',
				level => 2
			);

			@extraopts=();
			@extraopts=( rejectbade3 => 'N' )
			unless ($datamode eq "IMAGE.EL");

			# Check if an offset-variance map is available for this
			# instrument/CCD/node, if so add to the arguments for 'emevents'
			#
			my @ovfile=findFile(class => 'ODF',
					instrument => thisInstrument,
					ccd => $ccd,
					node => $node,
					content => 'Offset/variance mode events'
			);
			push(@extraopts, withoffvarsets => 'Y', offvarsets => join(" ",@ovfile))
			if @ovfile;
			ppd(id => "000006" , ccd => $ccd , node => $node) if @ovfile;

			doCommand('emevents',
				  odfeventset => $evFile1,
				  eventset => $evFile2,
				  analysepatterns => 'N',
				  flagbadpixels => 'N',
				  splitdiagonals => 'N',
				  withframeset => 'Y',
				  frameset => $frameFile,
				  randomizeposition => 'N',
				  @extraopts)
			or return exception();

			doCommand('emenergy',
				  ineventset => $evFile2,
				  correctcti => 'N',
				  correctgain => 'N',
				  randomizeenergy => 'N')
			or return exception();

			$fovArea += readFITSKeyword(file => $evFile2
						,extension => 'EVENTS'
						,keyword => 'IN_FOV'
			) if ( hasFITSKeyword(file => $evFile2
				,extension => 'EVENTS'
				,keyword => 'XMMEA_22'
				) 
			);
		}
    }

    # Make flare gti for use with bad pixel searching

    my $hasBadpixGti = 0;
    my $bpFlareGti;
    my $backgrndBinsize = 26;
    my $GTI_time_threshold=gtiTimeThreshold();

	my $fovAreaThreshold = 10;
    if( $fovArea > $fovAreaThreshold )
	{
		ppd(id => "000008" , data => { fovarea => $fovArea , threshold => $fovAreaThreshold } );
		my $bpEvFile=newFile(class => 'intermediate',
					 instrument => thisInstrument,
					 exp_id => $exp_id,
					 content => 'Bad pix find event list'
		);

		my @outerEvFiles=findFile(class => 'intermediate'
			,instrument => thisInstrument
			,exp_id => $exp_id
			,content => 'Processed events'
			,level => 2
		);

		doCommand('evlistcomb'
			,eventsets => [@outerEvFiles]
			,imagingset => $bpEvFile
			,instrument => 'emos'
			,emosdatamodes => 'IMAGING'
		);
		
		$bpFlareGti=newFile(class => 'intermediate',
					instrument => thisInstrument,
					exp_id => $exp_id,
					content => 'Bad pix flare GTI');
			my $histoFile=newFile(class => 'intermediate',
								  instrument => thisInstrument,
								  exp_id => $exp_id,
								  content => 'Flare rates file',
								  level => 1);

		if( hasFITSKeyword(file => $bpEvFile, extension => 'EVENTS', keyword => 'XMMEA_22') )
		{
			my $flareBitmask = getEventBitmask(file=>$bpEvFile, masks=>$flareFlags);
			ppd(id => "000009" , data => { flarebitmask => $flareBitmask } );

			doCommand('evselect'
				  ,table => $bpEvFile
				  ,writedss => 'N'
				  ,updateexposure => 'N'
				  ,withrateset => 'Y'
				  ,rateset => $histoFile
				  ,maketimecolumn => 'Y'
				  ,makeratecolumn => 'Y'
				  ,timebinsize => $backgrndBinsize
				  ,timecolumn => 'TIME'
				  ,expression => "(PATTERN==0) && #XMMEA_22 && ((FLAG & $flareBitmask) == 0)"
				  ,withimageset => 'N'
				  ) or return exception();
		
			if( numberFITSRows(file=>$histoFile, extension=>'RATE') > 0 && check_min_max_rate_col($histoFile) )
			{
# 				doCommand('tabcalc',
# 					  tables => "$histoFile:RATE",
# 					  columntype => 'real64',
# 					  column => 'RATE',
# 					  expression => "1000*COUNTS/($backgrndBinsize*$fovArea)" )
# 				or return exception();

# Apply relative corrections to flare rate file

# 				my $tmpHistoFile = newFile(class => 'intermediate'
# 							   , instrument => thisInstrument
# 							   , exp_id => $exp_id
# 							   , content => 'Copy flare rates file'
# 							   , level => 1);

# 				doCommand(
# 					  'epiclccorr'
# 					  , srctslist => $histoFile
# 					  , eventlist => $bpEvFile
# 					  , outset => $tmpHistoFile
# 					  , allcamera => 'Y'
# 					  , applyabsolutecorrections => 'N'
# 					  ) or return exception();

# 				copyFile(
# 					 source => $tmpHistoFile
# 					 , destination => $histoFile
# 					 );

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


				my $optCmd = "bkgoptrate fracexpstyle=auto tssettabname=$histoFile:RATE withendtime=N withlowercutoffcount=N withmintimeratio=N ignorelowcnttail=Y withstarttime=N -w 10";
				info("RUN: $optCmd");
				my $optRate = `$optCmd`;
				chomp $optRate;
				info("OPTRATE: $optRate");
				
				if ( $optRate ) {		# If $optRate is NOT empty

					doCommand('tabgtigen'
						  ,table => $histoFile
						  ,gtiset => $bpFlareGti
						  ,expression => "RATE<$optRate"
	#					  ,expression => 'RATE<2.0'
						  ,mingtisize => $mingtisize
					) or return exception();

					my $GTI_time=0;
					if ($GTI_time_threshold)
					{
						$GTI_time = readFITSKeyword(file => $bpFlareGti
									,extension => "STDGTI"
									,keyword => "ONTIME"
						);
						ppd(id => "000010" , data => { ontime => $GTI_time , threshold => $GTI_time_threshold } );
						info("Flare screening GTI_Time: $GTI_time, GTI_time_threshold: $GTI_time_threshold");
					}

					if  ( numberFITSRows(file =>  $bpFlareGti, extension => 'STDGTI') > 0
						&& ($GTI_time > $GTI_time_threshold || !$GTI_time_threshold)
					)
					{
						ppd( id => "000011" );
						$hasBadpixGti = 1;
					}
				} else {
					info("OPTRATE: Empty. Skipping 'Bad pix flare GTI' creation");
				}

				
			}
		}
    }

    $fovArea = 0;

    # Second loop over CCDs

    foreach my $evFile0 (@evFiles2) 
	{

		# Find CCD and node numbers
		my $ccd=readFITSKeyword(file => $evFile0,
					extension => 2,
					keyword => "CCDID");

                # Again, potentially skip some CCDs
                if (manualintervention( instrument => thisInstrument . "ccd", exposure => $ccd ))
                {
                    info("Manual Intervention. Ignoring CCD $ccd");
                    next;
                }

		my $node=readFITSKeyword(file => $evFile0,
					 extension => 2,
					 keyword => "CCDNODE");
		my $datamode=readFITSKeyword(file => $evFile0,
						 extension => 2,
						 keyword => "DATATYPE");
		info("Second pass of events for CCD $ccd, node $node, mode=$datamode");

		my $frameFile=findFile(class => 'intermediate',
					   instrument => thisInstrument,
					   exp_id => $exp_id,
					   node => $node, 'ccd' => $ccd,
					   content => 'Rectified frames');

		my $evFile1=findFile(class => 'intermediate',
					 instrument => thisInstrument,
					 exp_id => $exp_id,
					 node => $node,
					 ccd => $ccd,
					 content => 'Processed events',
					 level => 1);

		my $gtiFile1=findFile(class => 'intermediate',
					  instrument => thisInstrument,
					  exp_id => $exp_id,
					  node => $node,
					  ccd => $ccd,
					  content => 'CCD GTI',
					  level => 1);

		my $evFile3=newFile(class => 'intermediate',
					instrument => thisInstrument,
					exp_id => $exp_id,
					node => $node,
					ccd => $ccd,
					content => 'Processed events',
					level => 3);

		# Do bad pixel finding if imaging mode

		@extraopts=();

		if( $datamode eq "IMAGE.EL" || $datamode eq "RIMAGE.EL" )
		{
			my $evFile2=findFile(class => 'intermediate',
					instrument => thisInstrument,
					exp_id => $exp_id,
					node => $node,
					ccd => $ccd,
					content => 'Processed events',
					level => 2);

			my $badFile=newFile(class => 'intermediate',
					instrument => thisInstrument,
					exp_id => $exp_id,
					node => $node,
					ccd => $ccd,
					content => 'Bad pixel list');

			my $projFile=newFile(class => 'intermediate',
					instrument => thisInstrument,
					exp_id => $exp_id,
					node => $node,
					ccd => $ccd,
					content => 'Projected image');

			# Search for dark pixels, and bright pixels if no flare gti

			@extraopts=(findbright => 'N') if $hasBadpixGti;

			doCommand('emeventsproj',
				  eventset => $evFile2,
				  rejectbadevents => 'Y',
				  evimageset => $projFile)
			or return exception();

			doCommand('embadpixfind',
				  evimageset => $projFile,
				  badpixset => $badFile,
				  includedeadpixels => 'Y',
					  @extraopts)
			or return exception();

			# Search for bright pixels after filtering out flares if the falre gti is above the threshold

			if($hasBadpixGti)
			{
			
				doCommand('evselect'
					  ,writedss => 'N'
					  ,updateexposure => 'N'
					  ,table => $evFile2
					  ,keepfilteroutput => 'Y'
					  ,expression => "gti($bpFlareGti:STDGTI,TIME)"
					  ,withimageset => 'N'
				) or return exception();

				doCommand('emeventsproj',
					  eventset => $evFile2,
					  rejectbadevents => 'Y',
					  evimageset => $projFile)
					or return exception();

				doCommand('embadpixfind',
					  evimageset => $projFile,
					  badpixset => $badFile,
					  finddead => 'N',
					  incremental => 'Y')
					or return exception();
			}

			@extraopts=(getnewbadpix => 'Y'
				,badpixset => $badFile
				,getuplnkbadpix => 'Y'
				,getotherbadpix => 'Y'
			);
		}

		# run badpix on the event file
		#

		doCommand('badpix',
			  eventset => $evFile1,
			  withoutset => 'Y',
			  outset => $evFile3,
			  windowfilter => 'Y',
			  @extraopts)
			or return exception();

		# Check if an offset-variance map is available for this
		# instrument/CCD/node, if so add to the arguments for 'emevents'
		#
		@extraopts=();
		my @ovfile=findFile(class => 'ODF',
					instrument => thisInstrument,
					ccd => $ccd,
					node => $node,
					content => 'Offset/variance mode events');
	    if ( @ovfile )
		{
			ppd( id => "000012" );
			push(@extraopts, withoffvarsets => 'Y', offvarsets => join(" ",@ovfile));
		}


		#
		push(@extraopts, splitdiagonals => 'N')
			  if $datamode eq "TIME.EL" || $datamode eq "CTIME.EL";

		push(@extraopts, flagtruncatede1 => 'N')
			  if $datamode eq "CTIME.EL";

		#
		push(@extraopts, rejectbade3 => 'N')
		  unless $datamode eq "IMAGE.EL";


		my $evFile4=newFile(class => 'intermediate',
					instrument => thisInstrument,
					exp_id => $exp_id,
					node => $node,
					ccd => $ccd,
					content => 'Processed events',
					level => 4
		);

		doCommand("emevents",
			  withframeset => 'Y',
			  frameset => $frameFile,
			  odfeventset => $evFile3,
			  eventset => $evFile4,
			  randomizeposition => 'Y',
			  randomizetime => 'Y',
			  @extraopts
		) or return exception();


		# Align global gti file with event list and merge with ccd gti
			# Some files may not have been created, so need to check
			#
		my $hkGtiFileAligned=newFile(class => 'intermediate',
						 instrument => thisInstrument,
						 exp_id => $exp_id,
						 node => $node,
						 ccd => $ccd,
						 content => 'Aligned HK GTI');
		my $gtiFile2=newFile(class => 'intermediate',
					 instrument => thisInstrument,
					 exp_id => $exp_id,
					 node => $node,
					 ccd => $ccd,
					 content => 'CCD GTI',
					 level => 2);

		if( fileExists(file => $hkGTIFile) )
		{
			ppd("000013");
			if( ($datamode eq "IMAGE.EL" || $datamode eq "RIMAGE.EL") )
			{
				doCommand("gtialign",
					gtitable => $hkGTIFile . ":STDGTI",
				  	eventset => $evFile4,
				  	outset => $hkGtiFileAligned)
				or return exception();
			}
			else
			{
				copyFile(source => $hkGTIFile, destination => $hkGtiFileAligned); 
				doCommand("dscp",
					  from => "$hkGtiFileAligned:STDGTI",
					  to => "$hkGtiFileAligned:STDGTI" . $node . $ccd)
				  or return exception();
				doCommand("dsrm",
					  objects => "$hkGtiFileAligned:STDGTI")
				  or return exception();
			}

			doCommand("gtimerge",
				  tables => "$hkGtiFileAligned $gtiFile1",
				  mergemode => 'and',
				  gtitable => "$gtiFile2:STDGTI" . $node . $ccd)
			or return exception();
		}
		else
		{
			copyFile(source => $gtiFile1, destination => $gtiFile2);
		}

		# Add GTI file to evselect expression
		push(@filter,"(CCDNR==${node}${ccd} && GTI($gtiFile2,TIME))");

		# Apply attitude correction to events
		#
		doCommand("attcalc",
			  eventset => $evFile4,
			  refpointlabel => 'pnt',
			  withmedianpnt => 'true',
			  attitudelabel => 'ahf',
			  atthkset => $attFile)
			or return exception();


		# Rectify event energies
		#
		@extraopts=();

		# Background subtraction disabled in non-imaging mode
		push(@extraopts, getccdbkg => 'N', rejectbade3e4 => 'N')
		  unless $datamode eq "IMAGE.EL";

		doCommand('emenergy',
			  ineventset => $evFile4,
			  randomizeenergy => 'Y',
			  @extraopts)
		  or return exception();
		  
		$fovArea += readFITSKeyword(file => $evFile4
					,extension => 'EVENTS'
					,keyword => 'IN_FOV'
		) if ( hasFITSKeyword(file => $evFile4
			,extension => 'EVENTS'
			,keyword => 'XMMEA_22'
			) 
		);

		# Filter out unwanted events before lists are merged
		#
		my $eventBitmask = getEventBitmask(file=>$evFile4, masks=>$eventlistFlags);

		doCommand('evselect'
			  ,table => "${evFile4}:EVENTS"
			  ,destruct => 'Y'
			  ,keepfilteroutput => 'Y'
			  ,writedss => 'Y'
			  ,updateexposure => 'N'
			  ,expression => "(FLAG & $eventBitmask)==0"
			  ,withimageset => 'N'
		) or return exception();
    }


    # Merge per-CCD event lists into one per camera (for imaging and
    # timing mode)
    #
    @evFiles=findFile(class => 'intermediate',
		      instrument => thisInstrument,
		      exp_id => $exp_id,
		      content => 'Processed events',
		      level => 4);
    my $tmpimgEvents=newFile(class => 'intermediate',
			     instrument => thisInstrument,
			     exp_id => $exp_id,
			     content => 'imaging mode event list');
    my $tmptimEvents=newFile(class => 'intermediate',
			     instrument => thisInstrument,
			     exp_id => $exp_id,
			     content => 'timing mode event list');

    doCommand('evlistcomb', 
	      eventsets => [@evFiles],
	      imagingset => $tmpimgEvents, 
	      timingset => $tmptimEvents,
	      instrument => 'emos',
	      maintable => 'EVENTS OFFSETS',
	      othertables => 'EXPOSURE STDGTI BADPIX' )
	or return exception();


    # Create filter expression for event lists
    my $expr=join("||", @filter);


    # Find calibration index file
    my $cifFile=findFile(class => 'product',
			 content => 'Calibration index file');


    # Filter good imaging events into final product files
    if(fileExists(file => $tmpimgEvents))
	{
		my $imgEvents=newFile(class => 'product',
					  instrument => thisInstrument,
					  exp_id => $exp_id,
					  content => 'EPIC MOS imaging mode event list');
		my $eventBitmask = getEventBitmask(file=>$tmpimgEvents, masks=>$eventlistFlags);

		doCommand('evselect'
			  ,table => $tmpimgEvents
			  ,filteredset => "${imgEvents}:EVENTS"
			  ,withfilteredset => 'true'
			  ,expression => "($expr) && (FLAG & $eventBitmask)==0"
			  ,writedss => 'true'
			  ,cleandss => 'true'
			  ,updateexposure => 'true'
			  ,keepfilteroutput => 'true'
			  ,destruct => 'true'
			  ,withimageset => 'N'
		  ) or return exception();
    
    # Added to make the full image of a mosaic observation.
    
		doCommand('attcalc',
			  ,eventset      => $imgEvents,
		 	  ,refpointlabel => 'pnt',
	 		  ,withmedianpnt => 'true',
			  ,attitudelabel => 'ahf',
			  ,atthkset      => $attFile,
			  ,calctlmax     => 'yes')
			or return exception();

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

	# Find attitude GTI file
	    my $attGTIFile = findFile(
				      class => 'intermediate'
				      , instrument => 'all'
				      , content => 'attitude GTI'
				      );
	doCommand(
		'evselect'
		, table => $imgEvents
		, 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();

	# 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($imgEvents, "EVENTS", "TLMIN6" => {value => $tlmin6_inter});
	FITSUtils::writeKeyword($imgEvents, "EVENTS", "TLMAX6" => {value => $tlmax6_inter});
	FITSUtils::writeKeyword($imgEvents, "EVENTS", "TLMIN7" => {value => $tlmin7_inter});
	FITSUtils::writeKeyword($imgEvents, "EVENTS", "TLMAX7" => {value => $tlmax7_inter});
	#
	# 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
        #
        my $submode = readFITSKeyword(
				      file => $imgEvents
				      , extension => 'PRIMARY'
				      , keyword => 'SUBMODE'
				      ) or return exception();

        #
        # MOS case : always search for FWC data because peripheral CCDs are always in PrimeFullWindow mode
	#
	info("Imaging mode for ", thisInstrument , ", exposure: $exp_id : $submode. Producing FWC data...");

	#if ( $submode eq "PrimeFullWindow" ) {
	        #info("Imaging mode is PrimeFullWindow. Producing FWC data...");

                my $imgFWCEvents = newFile(
                    class => 'intermediate'
                    , instrument => thisInstrument
                    , exp_id => $exp_id
                    , content => 'EPIC MOS imaging mode FWC event list'
                );
		doCommand(
		    'evqpb'
		    , table => $imgEvents
		    , attfile => $attFile
		    , outset => $imgFWCEvents
		    , 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: $imgFWCEvents
		FITSUtils::writeKeyword($imgFWCEvents, "EVENTS", "TLMIN6" => {value => $tlmin6_inter});
		FITSUtils::writeKeyword($imgFWCEvents, "EVENTS", "TLMAX6" => {value => $tlmax6_inter});
		FITSUtils::writeKeyword($imgFWCEvents, "EVENTS", "TLMIN7" => {value => $tlmin7_inter});
		FITSUtils::writeKeyword($imgFWCEvents, "EVENTS", "TLMAX7" => {value => $tlmax7_inter});

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

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


		# Copy CIF into event list
		doCommand('dscopyblock',
			  blocks => "${cifFile}:CALINDEX",
			  to => $imgEvents)
		  or return exception();

		# Define flare gti
		# Works as soon as the GATTI was correctly reconstructed for one CCD
		#
		# only generate flareGTI given sufficient $fovArea

		if( $fovArea > 10 ){

		    # define flare GTI
		    my $flareGti=newFile(class => 'intermediate',
					 instrument => thisInstrument,
					 exp_id => $exp_id,
					 content => 'flare GTI');

		    # define flare BG timeseries
		    my $histoFile=newFile(class => 'intermediate',
					  instrument => thisInstrument,
					  exp_id => $exp_id,
					  content => 'Flare background timeseries');
		    
		    # check for XMMEA_22 keyword
		    if( hasFITSKeyword(file => $imgEvents, extension => 'EVENTS', keyword => 'XMMEA_22') ){
			
			my $flareBitmask = getEventBitmask(file=>$imgEvents, masks=>$flareFlags);

			# generate flare bg timeseries
			doCommand('evselect'
				  ,table => $imgEvents
				  ,writedss => 'Y'
				  ,updateexposure => 'N'
				  ,withrateset => 'Y'
				  ,rateset => $histoFile
				  ,timebinsize => $backgrndBinsize
				  ,maketimecolumn => 'Y'
				  ,timecolumn => 'TIME'
				  ,expression => "(PATTERN==0) && #XMMEA_22 && ((FLAG & $flareBitmask) == 0)"
				  ,withimageset => 'N'
				  ) or return exception();
			
			doCommand(
				  'tabcalc'
				  , tables => "$histoFile:RATE"
				  , columntype => 'real64'
				  , column => 'RATE'
				  , expression => "1000*COUNTS/($backgrndBinsize*$fovArea)"
				  )
			    or return exception();
			
			# Now subtract TSTART from the TIME value, leave the result in a new column T_ELAPSED:
			my $start_time=readFITSKeyword(
						       file => $histoFile
						       , extension => 'RATE'
						       , keyword => "TSTART"
						       );
			doCommand(
				  'tabcalc'
				  , tables => "$histoFile: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();

                        # Construct TS plot command file
			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', '')]
				       );
			
                        # Generate flare background PS
			my $flareBkgPs = newFile(
						 class => 'intermediate'
						 , instrument => thisInstrument
						 , exp_id => $exp_id
						 , content => 'F bkg ts'
						 , format => 'PS'
						 );
			doCommand(
				  'fplot'
				  , infile => $histoFile
				  , xparm => 'T_ELAPSED'
				  , yparm => 'RATE'
				  , rows => '-'
				  , device => "$flareBkgPs/cps"
				  , pltcmd => '@'."$fplotCommandFile"
				  )
			    or return exception();

			my $mos_tbinsize;
			my $sourceFreeBackRates;
			my $SFBLCRegionFileFaint;
			my $SFBLCRegionFileBright;
			my $SFBLCBackScal;

			if( numberFITSRows(file=>$histoFile, extension=>'RATE') > 0 && check_min_max_rate_col($histoFile) )
			{

# run bkgoptrate for first determination of optimal flare cutoff			    
			    my $retval = 0;
			    my $optCmd = "bkgoptrate fracexpstyle=auto tssettabname=$histoFile:RATE withendtime=N withlowercutoffcount=N withmintimeratio=N ignorelowcnttail=Y withstarttime=N -w 10";
			    info("CMD: $optCmd");
			    my $optRate = `$optCmd`;
			    $retval = $?;
			    if ($retval ne 0) { return exception("Error $retval returned by bkgoptrate: $!"); }

			    chomp $optRate;
			    info("OPTRATE: $optRate");
			    
			    info("Performing second optimised GTI filter < $optRate");
			    doCommand('tabgtigen'
				      ,table => $histoFile
				      ,gtiset => $flareGti
				      ,expression => "RATE<$optRate"
				      ,mingtisize => $mingtisize
				      ) or return exception();
			    
			}

			info("==== starting source-free background flare lightcurve processing===");
			
                        # Generate image for SFBLC analysis

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

			# Find attitude GTI file
			my $attGTIFile = findFile(
						  class => 'intermediate'
						  , instrument => 'all'
						  , content => 'attitude GTI'
						  );

                        # define MOS filter expression
			my $expr;
			if (fileExists( file => $flareGti ) ){
			    $expr = "(GTI($flareGti, TIME) && GTI($attGTIFile, TIME) && PATTERN <= 12 && PI in (500:7500] && (FLAG & 0x766aa000)==0)";
			} else {
			    $expr = "(GTI($attGTIFile, TIME) && PATTERN <= 12 && PI in (500:7500] && (FLAG & 0x766aa000)==0)";
			}

                        # makeimage (make use of existing flare GTI)
			# In DET coordinates

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

			doCommand(
				  'evselect'
				  , table => $imgEvents
				  , imageset => $SFBLCInputImage
				  , xcolumn => 'DETX'
				  , ycolumn => 'DETY'
				  , imagebinning => 'binSize'
				  , ximagebinsize => $xbin_SFBLC
				  , yimagebinsize => $ybin_SFBLC
				  , expression => $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_ts_expr;
#			if (! &ModuleUtil::isImageLimit($SFBLCInputImage,'',1000.0)) {
#
#			The Mask is produced although image counts < 1000
#
			    # make vignetted exposure map
			    
			    info("Creating 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 => $imgEvents
				      , 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 => $imgEvents
				      , attitudeset => $attFile
				      , attrebin => 2
				      , expimageset => $SFBLCNonVigExpmap
				      , pimin => 500
				      , pimax => 7500
				      , withvignetting => 'no'
                                      , withdetcoords => 'yes'
				      ) 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
				      ) or return exception();
			    
##########################################################################################################################################
			# check at least 1000 counts in input image before continuing to process: Source detection, etc

			my $SFBLCDetMaskArith;

			my $source_free_ts_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 MOS, 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'
						       );
			    
			    my $SFBLCBoxListFaint = newFile(
							    class => 'intermediate'
							    , instrument => thisInstrument
							    , exp_id => $exp_id
							    , content => 'SFBLC box-local source list faint subset'
							    , band => 8
							    , format => 'FITS'
							    );
			    
			    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();
			    

			    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();
			    }
			    
			    # Now need to isolate ID_INT=ID_BAND=0 lines to a new file
			    # and convert to a region file
			    
			    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'
							      );
			    }
			    
			    # Generate region selection expression
			    
			    my $lthresh = 50; # ML threshold for removing bright sources  
			    my $rad = 60; # Fixed extraction radius in arcsec (faint)
			    my $rad2 = 100; # Fixed extraction radius in arcsec (bright)
			    my $region_select;
			    my $region_select_drop;
			    if (thisInstrument =~ /mos1/) {
				$region_select = "(ID_BAND == 1 && ID_INST == 2 && LIKE >= $lthresh && BOX_CTS > 0)";
			    } else {
				$region_select = "(ID_BAND == 1 && ID_INST == 3 && LIKE >= $lthresh && BOX_CTS > 0)";
			    };
			    
			    # For source masking (drop)
			    if (thisInstrument =~ /mos1/) {
				$region_select_drop = "(ID_BAND == 1 && ID_INST == 2 && LIKE >= 30 && BOX_CTS > 0)";
			    } else {
				$region_select_drop = "(ID_BAND == 1 && ID_INST == 3 && LIKE >= 30 && BOX_CTS > 0)";
			    };
			    
			    doCommand(
				      'region'
				      , eventset => "${imgEvents}:EVENTS"
				      , tempset => $SFBLCTempFile
				      , srclisttab => $SFBLCBoxListFaint
				      , expression => $region_select
				      , operationstyle => 'global'
				      , fixedradius => $rad
				      , bkgregionset => "${SFBLCRegionFileFaint}:REGION"
				      , outunit => "detxy"
				      ) or return exception();
			    
			    doCommand(
				      'region'
				      , eventset => "${imgEvents}: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 => "${imgEvents}: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 $sourceFreeEvalImage = newFile(
							      class => 'intermediate'
							      , instrument => thisInstrument
							      , exp_id => $exp_id
							      , content => 'SFBLC evaluation image'
							      , band => 8
							      , format => 'FITS'
							      );
			    
			    # include GTI, regions in source-free image selection expression

			    my $source_free_img_expr;
			    if ( fileExists( file => $flareGti ) ){
			        $source_free_img_expr = "((GTI($flareGti, TIME) && GTI($attGTIFile, TIME)) && PATTERN<=12 && PI in (500:7500] && (FLAG & 0x766aa000)==0 && (region($SFBLCRegionFileFaint,DETX,DETY) && region($SFBLCRegionFileBright,DETX,DETY)))";
			    } else {
			        $source_free_img_expr = "((GTI($attGTIFile, TIME)) && PATTERN<=12 && PI in (500:7500] && (FLAG & 0x766aa000)==0 && (region($SFBLCRegionFileFaint,DETX,DETY) && region($SFBLCRegionFileBright,DETX,DETY)))";
			    }

			    # In DET coordinates
			    doCommand(
				      'evselect'
				      , table => $imgEvents
				      , imageset => $sourceFreeEvalImage
				      , xcolumn => 'DETX'
				      , ycolumn => 'DETY'
				      , withimageset => 'Y'
				      , imagebinning => 'binSize'
				      , ximagebinsize => $xbin_SFBLC
				      , yimagebinsize => $ybin_SFBLC
				      , expression => $source_free_img_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
			    
			    # MOS
			    # SRR thinks input event file is already modified to include CCD GTIs and
			    # badpixel info. So here can just run evselect, expression is modified
			    # to exclude XMMEA_22 and includes the energy range
			    
			    # include GTI, regions in source-free time-series expression

			    # try merging region subsets as input to evselect

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

#			    $source_free_ts_expr = "PATTERN<=12 && PI in (500:7500] && ((FLAG & 0x762ba000)==0) && (region($SFBLCRegionFileFaint,X,Y) && region($SFBLCRegionFileBright,X,Y))";
			    $source_free_ts_expr = "PATTERN<=12 && PI in (500:7500] && ((FLAG & 0x762ba000)==0) && (region($SFBLCRegionFileMerged,DETX,DETY))";

			} else {

			    # extraction expression for sparse images (omit region select)
			    info("Creating in-band light curve for MOS (sparse image < 1000 counts)");

			    $source_free_ts_expr = "PATTERN<=12 && PI in (500:7500] && ((FLAG & 0x762ba000)==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]";

			};

		#
		#  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'
						);
		#info("DEBUG - xbin, ybin = $xbin , $ybin");

		$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 time series as product

			$mos_tbinsize = 26;
			$sourceFreeBackRates = newFile(
						       class => 'product'
						       , instrument => thisInstrument
						       , exp_id => $exp_id
						       , content => 'EPIC flare background timeseries'
						       );
			
			doCommand(
				  'evselect'
				  , table => $imgEvents
				  , expression => $source_free_ts_expr
				  , maketimecolumn => 'Y'
				  , makeratecolumn => 'Y'
				  , rateset => $sourceFreeBackRates
				  , timebinsize => $mos_tbinsize
				  , timecolumn => 'TIME'
				  , updateexposure => 'N'
				  , withrateset => 'Y'
				  , writedss => 'Y'
				  ) or return exception();

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

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


			# Check out EXPOSURE keyword in $sourceFreeBackRates
			#
			my $optFlareGti;
			my $sourceFreeLCCommandFile;
			if ( readFITSKeyword( file => $sourceFreeBackRates, extension => 'RATE', keyword => 'EXPOSURE' ) ){
			
				doCommand(
					  'epiclccorr'
					  , srctslist => $sourceFreeBackRates
					  , eventlist => $imgEvents
					  , 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 $hasOptFlareGti = 0;
					my $rcut_mos;

					if( numberFITSRows(file=>$sourceFreeBackRates, extension=>'RATE') > 0 && check_min_max_rate_col($sourceFreeBackRates) ){

						my $retval = 0;
						$rcut_mos = `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: $!"); }

						if ($rcut_mos eq '') {
						    info("bkgoptrate call returns no value, skipping optimised flare GTI generation...");
						} else {

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

							info("bkgoptrate returns optimum cut threshold: $rcut_mos");

							# append SN_TO_BKGRATE 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_mos
								  ) or return exception();

							# clamp rcut to maximum 100
							#
							#if ($rcut_mos > 100) {
							#	$rcut_mos = 100;
							#	info("optimum cut rate threshold reset 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_mos}"
								  , mingtisize => 100
								  ) or return exception();

						    }
						$hasOptFlareGti = 1;
						}
					} else {
					   info("Flat TS. Do not run bkgoptrate, skipping optimised flare GTI generation...");
					}



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

					$sourceFreeLCCommandFile = newFile(
									      class => 'intermediate'
									      , instrument => thisInstrument
									      , exp_id => $exp_id
									      , content => 'In-band sf lc pco file'
									      , format => 'ASCII'
						      , extension => 'pco'
						      );

					# Cutoff threshhold label string
					my $srcut = sprintf("%.3f",$rcut_mos);

					# 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 ($snMin, $snMax);
					if (fileExists( file => $snDiagFile )){
					my $snBg = readFITSColumn(
								  file => $snDiagFile
								  , extension => 'SN_TO_BKGCUT'
								  , column => 'SN_RATIO'
								  ) or return exception();

					$snMin = min(@$snBg);
					$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 noice 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();

					doCommand(
						  'fplot'
						  , infile => "${sourceFreeBackRates}[SN_TO_BKGCUT]"
						  , 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'
							      );

					my $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();


			} else {
				info("No EXPOSURE keyword in TS, then No 'EPIC flare background timeseries'");
			}

			
                        # 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 = "GTI($attGTIFile, TIME) && PATTERN<=12 && PI in (500:7500] && ((FLAG & 0x762ba000)==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 => $imgEvents
				  , timebinsize => $mos_tbinsize
				  , updateexposure => 'N'
				  , withrateset => 'Y'
				  , writedss => 'Y'
				  ) or return exception();
			
# 			doCommand(
# 				  'tabcalc'
# 				  , column => 'RATE'
# 				  , columntype => 'real64'
# 				  , expression => "COUNTS/${mos_tbinsize}"
# 				  , tables => "${optBackRates}:RATE"
# 				  ) or return exception();
			
                        # Apply relative corrections to test GTI ts via epiclccorr()

			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( (fileExists(file => $optFlareGti) && numberFITSRows(file=>$optFlareGti, extension=>'STDGTI') > 0) &&
                            !ModuleUtil::isColumnEmpty($optBackRates,"RATE") )
			{
			    
			    doCommand(
				      'epiclccorr'
				      , srctslist => $optBackRates
				      , eventlist => $imgEvents
				      , 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();
			
			}
		    }
		}
	    }
	
	
	# Filter good timing events into final product files
	if(fileExists(file => $tmptimEvents)){

	    my $timEvents=newFile(class => 'product',
				  instrument => thisInstrument,
				  exp_id => $exp_id,
				  content => 'EPIC timing mode event list');

	    doCommand('evselect'
		      ,table => $tmptimEvents
		      ,filteredset => $timEvents
		      ,withfilteredset => 'true'
		      ,keepfilteroutput => 'true'
		      ,updateexposure => 'true'
		      ,writedss => 'true'
		      ,cleandss => 'true'
		      ,expression => $expr
		      ,destruct => 'true'
		      ,withimageset => 'N'
		      ) or return exception();
	    
	    # Copy CIF into event list
	    doCommand('dscopyblock',
		      blocks => "${cifFile}:CALINDEX",
		      to => $timEvents)
		or return exception();
	}
	
	return success();
    }

# -----------------------------------------------------------------------------
sub check_min_max_rate_col
{
    my ($file) = @_;

    my $fileRATE = readFITSColumn(
       file => $file
       , extension => 'RATE'
       , column => 'RATE'
       ) or return exception();

       my $fileRATEMin = min(@$fileRATE);
       my $fileRATEMax = max(@$fileRATE);

    if ( $fileRATEMin != $fileRATEMax ){
	return 1;
    } else { return 0; }

}
# -----------------------------------------------------------------------------

1;