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

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

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

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

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

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

# 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`;
			    chomp $optRate;
			    $retval = $?;
			    if ($retval ne 0) { return exception("Error $retval returned by bkgoptrate: $!"); }
			    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 = "(GTI($flareGti, TIME) && GTI($attGTIFile, TIME) && PATTERN <= 12 && PI in (500:7500] && (FLAG & 0x766aa000)==0)";

                        # makeimage (make use of existing flare GTI)

			doCommand(
				  'evselect'
				  , table => $imgEvents
				  , imageset => $SFBLCInputImage
				  , xcolumn => 'X'
				  , ycolumn => 'Y'
				  , imagebinning => 'binSize'
				  , ximagebinsize => 80
				  , yimagebinsize => 80
				  , 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)) {
			    
			    # 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'
				      ) 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'
				      ) 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();
			    
			    
			    # 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'
							     );

			    # 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();
			    
			    # 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'
							     );
			    
			    # 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;
			    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)";
			    };
			    
			    doCommand(
				      'region'
				      , eventset => "${imgEvents}:EVENTS"
				      , tempset => $SFBLCTempFile
				      , srclisttab => $SFBLCBoxListFaint
				      , expression => $region_select
				      , operationstyle => 'global'
				      , fixedradius => $rad
				      , bkgregionset => "${SFBLCRegionFileFaint}:REGION"
				      ) or return exception();
			    
			    doCommand(
				      'region'
				      , eventset => "${imgEvents}:EVENTS"
				      , tempset => $SFBLCTempFile
				      , srclisttab => $SFBLCBoxListBright
				      , expression => $region_select
				      , operationstyle => 'global'
				      , fixedradius => $rad2
				      , bkgregionset => "${SFBLCRegionFileBright}:REGION"
				      ) 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 = "((GTI($flareGti, TIME) && GTI($attGTIFile, TIME)) && PATTERN<=12 && PI in (500:7500] && (FLAG & 0x766aa000)==0 && (region($SFBLCRegionFileFaint,X,Y) && region($SFBLCRegionFileBright,X,Y)))";
			    
			    doCommand(
				      'evselect'
				      , table => $imgEvents
				      , imageset => $sourceFreeEvalImage
				      , xcolumn => 'X'
				      , ycolumn => 'Y'
				      , withimageset => 'Y'
				      , imagebinning => 'binSize'
				      , ximagebinsize => 80
				      , yimagebinsize => 80
				      , 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,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();


#			    $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,X,Y))";
			    
			} 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)";
			    			    
			};
			
			# 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'
							     );

			doCommand(
				  'epiclccorr'
				  , srctslist => $sourceFreeBackRates
				  , eventlist => $imgEvents
				  , 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_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: $!"); }

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

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

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

1;