package pnEvents;
use strict;
use English;
use Carp;
use List::Util qw(max min);
use Regexp::Common qw(boolean);
use vars
qw(@ISA $VERSION $name $author $date $version @instList $numberOfStreams);
@ISA = qw(Module);
use ModuleResources;
# Declare identity, version, author, date, etc.
$name = __PACKAGE__;
$VERSION = '2.46';
$version = $VERSION;
$author = 'Dean Hinshaw,Richard West,Duncan John Fyfe,Ian Stewart,Duncan Law-Green,Ed Chapin';
$date = '2013-08-26';
#
# ChangeLog
# =========
#
# Version 2.46 - 2013-08-26 (EC)
# ------------
#
# + RATE column from evselect must be non-zero for final epiclccorr call
#
# Version 2.45 - 2013-06-11 (EC)
# ------------
#
# + Merge changes from SOC version.
#
# version 2.44 - 2013-01-03 DLG
#
# - shorten intermediate SFBLC region file subset names to match MOS convention
#
# version 2.43 - 2012-12-03 DLG
#
# - merge faint and bright region files for generation of source-free background light curve
#
# version 2.42 - 2012-11-06 DLG
#
# - change references to $pn_timebinsize to $timebinsize
#
# version 2.41 - 2012-10-25 DLG
#
# - require valid nonzero flare GTI for final epiclcorr call
#
# version 2.40 - 2012-10-25 DLG
#
# - Modify plotting of SN output graph
#
# version 2.39 - 2012-10-24 DLG
#
# - evselect makeratecolumn=Y, disable tabcalc cal
#
# version 2.38 - 2012-10-23 DLG
#
# - added allcamera=Y to epiclccorr calls
#
# version 2.37 - 2012-10-17 DLG
#
# - check for existence of SFBLC region files in evselect
#
# version 2.36 - 2012-10-09 DLG
#
# - set writedss=Y on remaining SSP evselect calls
# = remove writedss param from epiclccorr
#
# version 2.35 - 2012-10-08 DLG
#
# - try writedss=Y on epiclccorr call
# + move tabcalc T_ELAPSED after epiclccorr
#
# version 2.34 - 2012-10-03 DLG
#
# - reenable epiclccorr, specify outset as temp file
#
# version 2.33 - 2012-10-01 DLG
#
# - disable epiclccorr calls pending GTI crash fix
#
# version 2.32 - 2012-09-27 DLG
#
# + correct epiclccorr param to 'applyabsolutecorrections'
#
# version 2.31 - 2012-09-26 DLG
#
# + Apply relative corrections to PN timeseries via epiclccorr()
#
# version 2.30 - 2012-09-26 DLG
#
# + undo 2.29 (moved to EPICSourceProducts)
# + correct evselect logic for region files
# + modify RATE column to real64
#
# version 2.29 - 2012-09-04 DLG
#
# + copy region file extensions to FBKTSR product
#
# version 2.28 - 2012-09-04 DLG
#
# + copy region file extensions to FBKTSR product
#
# version 2.27 - 2012-09-04 DLG
#
# + Add source regions to evselect for Optimised flare bkg time series
#
# version 2.26 - 2012-08-15 DLG
#
# + append SN_TO_BKGCUT table on background time series product
#
# version 2.25 - 2012-07-27 DLG
#
# + modify SFBLC region code: source exclusion radius to be dependent on source counts
#
# version 2.24 - 2012-07-23 DLG
#
# + disable 100cts limit on $rcut_pn
#
# version 2.23 - 2012-07-18 DLG
#
# + add call to epnoise to flag noisy CCDs
#
# version 2.22 - 2012-07-09 DLG
#
# + add ignorelowcnttail=yes parameter for bkgoptrate calls
#
# version 2.21 - 2012-07-02 DLG
#
# + Refactor code for sparse image handling
#
# version 2.20 - 2012-06-25 DLG
#
# + Remove global attitude GTI filter constraint from in-band SF light curve
#
# version 2.19 - 2012-06-25 DLG
#
# + Add FLCUTTHR header keyword to opt flare light curve
#
# version 2.18 - 2012-06-22 DLG
#
# + Do not perform source-free background analysis for sparse images < 1000 counts
#
# version 2.17 - 2012-06-21 DLG
#
# + Generate in-band source-free light curves FBKTSR as products (FITS, PDF)
#
# version 2.16 - 2012-06-20 DLG
#
# + Optimal cut threshold limited to 100 for flare GTI generation
#
# version 2.16 - 2012-05-25 DLG
#
# + Add SN diagnostic for source-free bkgoptrate calc
#
# version 2.15 - 2012-05-25 DLG
#
# + Add source-free flare background extraction
#
# version 2.14 - 2012-04-26 DLG
#
# + Merge epexposure call from saseval8.4
#
# version 2.13 - 2011-03-07 DLG
#
# + Add test call to bkgoptrate, start optimised flare background subtraction development
#
# version 2.12 - 2009-10-05 DLG
#
# + Plot EPIC background flare light curves in log space (Mantis 52)
#
# version 2.11 - 2009-09-20 DLG
#
# + Added counting cycle report check, set $mingtisize accordingly
#
#
# version 2.10 - 2009-07-28 DLG
#
# + set mingtisize=0 for testing
#
# version 2.09a - 2009-06-19 DLG
#
# + ignored if SAS_ATTITUDE set to "RAF" for slew processing
#
# version 2.09 - 2009-06-11 DLG
# ------------
#
# + epreject calls on imaging mode only per Mantis #46
#
#
# version 2.08 - 2009-05-15 DLG
# ------------
#
# + added support for epreject per Mantis feature req #35
#
#
# version 2.07 - 2009-01-12 DLG
# ------------
#
# + Changed "withtempcorrection" to "Y" following advice from Hermann Brunner
#
#
# version 2.06 - 2006-07-10 DJF
# ------------
#
# + Fix for darck quadrant problem applied to OAL. No longer need OBT_WARN test.
# Replaced test with an echo of the OBT_WARN value to the log file.
#
# version 2.05 - 2006-04-04 DJF
# ------------
#
# + Test for OBT_WARN keyword in epframes output. Return ignore if it is found to be true.
#
# version 2.04 - 2005-11-16 IMS
# ------------
#
# + Incorporated DJF's shortened intermediate file names.
#
# version 2.03 - 2005-10-14 IMS
# ------------
#
# + TSTART obtained from flare time series header and subtracted from the times of the flare PDF plot.
# + Flare plot Y axis now RATE not COUNTS. Calculation of RATE is thus moved forward a bit, and all flare productions now moved within the test for existence of RATE extension.
#
# version 2.02 - 2005-08-12 RGW
# ------------
# + changes to flare plotting to work around 80-character filename limit
# in 'fplot'
#
# version 2.01 - 2005-08-02 IMS
# ------------
#
# + PDF product now made from flare bkg time series.
# + Slightly changed event selection expression (events with PI < 150 eV are no longer kept; also the flag filter mask is given explicitly, not as XMMEA_EP).
#
# version 2.00 - 2005-05-09 DJF
# ------------
#
# + Add explicit perl module headers. Previously these were supplied
# at compile time. This will make debugging and extending the modules
# through additional perl libraries easier.
#
# + Code layout made more uniform with perltidy
# + Standerdized date format (ccyy-mm-dd)
# + Standerdized author list to use comma separator
# + Make use of perl $VERSION magic. Now $Version = version = ''
##
# Version 1.69 - 16-Mar-2004 (DJF)
# ------------
#
# + Added dummy badpixset=temporary file to first call to badpixfind
# to stop it creating it's default.
#
# Version 1.68 - 14-Jan-2004 (DJF)
# ------------
#
# + Use tabgtigen parameter mingtisize rather than evselect to exclude
# small GTI from flare background
#
# Version 1.67 - 10-12-2003 (DJF)
# ------------
#
# + Adapted doCommand to use anonymous lists for list parameters
# + Fixed filter expression. PI<-150 should have been PI<=150
#
# Version 1.66 - 05-Sept-2003 (DJF)
# -------------------------------
# + Fixed backgroundrate in second call to badpixfind. Was 0.0011 should have been 0.00011
#
# Version 1.65 - 11-Nov-2002 (DJF)
# -------------------------------
#
# + Updated evlistcomb columns to include OFFSETX
# + Added DLIMAP to evlistcomb othertables option
#
# Version 1.63 - 11-Nov-2002 (DJF)
# -------------------------------
#
# + Updated backgroundrate (0.0001 -> 0.0011) and hienergythresh (10.0 -> 12.0)
# in second call to badpixfind
#
# Version 1.62 - 17-Oct-2002 (DJF)
# -------------------------------
#
# + Moved OFFSETS from othertables to maintable option
#
# Version 1.61 - 19-Feb-2002 (DH)
# -------------------------------
#
# + Add OFFSETS to the list of tables for evlistcomb to propagate.
#
# Version 1.60 - 15-Feb-2002 (DH)
# -------------------------------
#
# + Bug fixes for 1.60 .
#
# Version 1.59 - 7-Feb-2002 (DH)
# -------------------------------
#
# + Remove RAWY>10 cut on event lists in imaging mode.
#
# Version 1.58 - 6-Nov-2001 (DH)
# -------------------------------
#
# + New parameters for badpixfind for making of the background
# mask. Values as recommended by Andy Read, via email from Michael
# Freyburg dated 10/10/01.
#
# Version 1.57 - 9-Oct-2001 (DH)
# -------------------------------
#
# + Change energy thresholds in badpixfind for making of the background
# mask. Was 0.16/10.0, now 7.0/15.0 .
#
# Version 1.56 - 4-Oct-2001 (DH)
# -------------------------------
#
# + Fix bug that made flare screening only work for full frame modes.
# + Change threshold on number for good background pixels from 0 to 100.
#
# Version 1.55 - 2-Oct-2001 (DH)
# -------------------------------
#
# + Put back flare gti creation for small window mode. Instead
# inhibit flare gti creation if the number of suitable background
# pixels found in a field is zero.
# + Increase flare gti cut from 6 to 10.
#
# Version 1.54 - 1-Oct-2001 (DH)
# -------------------------------
#
# + Scale integration time for flare histograms ($timebinsize)
# to 2 times nominal for large window mode, and 35 times for
# small window mode.
# + Do not create flare gti file for small window mode, though
# a flare histogram is still created.
#
# Version 1.53 - 1-Oct-2001 (DH)
# -------------------------------
#
# + Replace dslatts call with the hasFITSExtension function.
#
# Version 1.52 - 2001-08-09 (DJF)
# ------------
#
# + Changed Flare light curve from intermediate to final product.
#
# Version 1.51 - 2001-06-15
# ------------
#
# + Change withtempcorrection parameter back to 'N'.
#
# Version 1.50 - 2001-06-06
# ------------
#
# + Change withtempcorrection=N parameter of epevents to 'Y'.
#
# Version 1.49 - 2001-06-05
# ------------
#
# + Re-work merging of HK and CCD specific gti data. gtialign
# was not being called properly.
#
# Version 1.48 - 2001-05-22
# ------------
#
# + Only apply RAWY>10 only for imaging mode data.
# + Add withtempcorrection=N parameter to epevents, in case the
# default changes.
#
# Version 1.47 - 2001-05-16
# ------------
#
# + Bug fixes and refinements for changes made in version 1.46
#
# Version 1.46 - 2001-05-16 (DH)
# ------------
#
# + Remove unwanted columns from the EXPOSUnn extensions
# of the pn event lists.
#
# Version 1.45 - 2001-04-17 (DH)
# ------------
#
# + Changes to main badpixfind call:
# - Change loenergythresh from 0.16 to 0.14
# - Change hicolthresh from 0.0015 to 0.00105
#
# Version 1.44 - 2001-03-16 (DH)
# ------------
#
# + Print out version number in performAction() for
# tracking purposes.
#
# Version 1.43 - 2001-02-01 (DH)
# ------------
#
# + Change threshold for flare rejection rate cutoff from 36.4 to 6
#
# Version 1.42 - 2001-01-26 (DH)
# ------------
#
# + Restrict event energy range for flare rejection to 7-15 keV
#
# Version 1.41 - 2001-01-19 (DH)
# ------------
#
# + Change narrowerthanpsf parameter in standard call
# to badpixfind from 1.25 to 1.5.
#
# Version 1.40 - 2001-01-19 (DH)
# ------------
#
# + hicolthesh values for the two calls to badpixfind
# were not set correctly. In particular, value for the
# standard call was off by a factor of 10. Values now
# set as per Andy Read's email dated 13-Dec-2000.
#
# Version 1.39 - 2001-01-11 (DH)
# ------------
#
# + Remove GlobalHK from ignore rule.
#
# + Add InstrumentHK to start rule.
#
# Version 1.38 - 2000-12-19 (DH)
# ------------
#
# + Change hicolthesh in badpixfind masking back
# to its nominal value.
#
# Version 1.37 - 2000-12-13 (DH)
# ------------
#
# + First production version.
#
# Declare list of instruments this module is interested in
@instList = qw(epn);
# Number of streams
sub numberOfStreams
{
return numberOfExposures();
}
# Rules method
sub evaluateRules
{
#
return ignore()
if ignored(
module => 'ExpChooser'
, instrument => thisInstrument
, stream => thisStream
);
# ignore for slew processing
return ignore() if ($ENV{'PCMS_ISSLEW'});
#
start()
if complete(
module => 'ExpChooser'
, instrument => thisInstrument
, stream => thisStream
)
and complete(
module => 'InstrumentHK'
, instrument => thisInstrument
, stream => 1
)
and complete(
module => 'GlobalHK'
, instrument => 'all'
, stream => 1
);
}
# Action method
sub performAction
{
info("Module version number: $version");
my @extraopts;
# Work out which exposure this stream maps on to
my $exp_id = exposureID(
instrument => thisInstrument
, stream => thisStream
);
info("Exposure ID: $exp_id");
#
my $filter = getExposureProperty(
exp_id => $exp_id
, name => 'filter'
);
my $mode = getExposureProperty(
exp_id => $exp_id
, name => 'mode'
);
info("Filter: $filter");
info("Instrument mode: $mode");
# Find event files for this exposure
my @evFiles = findFile(
class => 'ODF'
, instrument => thisInstrument
, exp_id => $exp_id
, content => '*events'
);
# Return immediately if there are no event lists
return success() if $#evFiles < 0;
# Find attitude HK file
my $attFile = findFile(
class => 'product'
, instrument => 'all'
, content => 'Attitude time series'
, required => 'true'
);
# Find HK GTI file
my $hkGTIFile = findFile(
class => 'intermediate'
, instrument => thisInstrument
, content => 'HK GTI'
);
# Set mingtisize to 100 by default
my $mingtisize = 100;
# Check for counting cycle report
# reset mingtisize if required
my @ccxfiles = findFile(class => 'ODF'
,instrument => thisInstrument
,exp_id => $exp_id
,content => 'Counting cycle report'
);
unless (@ccxfiles) {
info("No counting cycle report found. Setting mingtisize=0");
$mingtisize = 0;
}
#
my ( @finalEvents, @gtiFilter, @maskFilter, $datamode );
my $goodPixels = 0;
# Time binning to use when making background light curves
#
my $timebinsize = 10.0;
# Loop through each event file
foreach my $evFile0 (@evFiles)
{
# Find CCD and node numbers
my $ccdid = readFITSKeyword(
file => $evFile0
, extension => 2
, keyword => "CCDID"
);
my $quadrant = readFITSKeyword(
file => $evFile0
, extension => 2
, keyword => "QUADRANT"
);
my $ccd = $quadrant * 3 + $ccdid + 1;
$datamode = readFITSKeyword(
file => $evFile0
, extension => 2
, keyword => "DATATYPE"
);
# my $ccdmode = readFITSKeyword(
# file => $evFile0
# , extension => 2
# , keyword => "CCDMODE"
# );
# info("Processing events for CCD $ccd, mode $datamode.");
# Unused. IMS 2005-08-03
my $evFile1 = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, ccd => $ccd
, content => 'Processed events'
, level => 1
);
my $gtiFile1 = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, ccd => $ccd
, content => 'CCD GTI'
, level => 1
);
# Run epframes
#
doCommand(
'epframes'
, set => $evFile0
, wrongpixlimit => 20
, srcposition => 190
, eventset => $evFile1
, gtiset => $gtiFile1
, mipmethod => 'onboard'
, mipthreshold => 3000
, setupbpx => 'nom3'
)
or return exception();
my %obt_warn = ( file => $evFile1
, extension => 'EXPOSURE'
, keyword => 'OBT_WARN'
);
if (hasFITSKeyword( %obt_warn ))
{
my $opt_warn = readFITSKeyword( %obt_warn );
if ( $opt_warn =~ /$RE{boolean}{true}/ )
{
info("OBT_WARN True");
}
else
{
info("OBT_WARN False");
}
}
# Run badpixfind to create a mask to be used for
# generating a background lightcurve
#
my $dummyBadPixels = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, ccd => $ccd
, content => 'Ignore This Bad pixels'
);
my $badPixMask = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, ccd => $ccd
, content => 'Bright pixels mask'
);
doCommand(
'badpixfind', eventset => $evFile1
, badpixmap => $badPixMask
, badpixset => $dummyBadPixels
, withfovmask => 'Y'
, threshabovebackground => 'Y'
, hithresh => 0.00025
, hicolthresh => 0.00007
, narrowerthanpsf => 0
, withbadpixmap => 'Y'
, withsearchbadcolumn => 'Y'
, loenergythresh => 7.0
, hienergythresh => 15.0
, backgroundrate => -1
, thresholdlabel => 'rate'
)
or return exception();
$goodPixels += readFITSKeyword(
file => $badPixMask
, extension => 1
, keyword => "GOODPIX"
);
push( @maskFilter
, "(CCDNR==${ccd} && MASK($badPixMask,0,0,RAWX,RAWY))" );
# Run badpixfind
#
@extraopts = ();
if ( $datamode eq "IMAGE.EL" )
{
my $newBadPixels = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, ccd => $ccd
, content => 'Bad pixels'
);
doCommand(
'badpixfind'
, eventset => $evFile1
, badpixset => $newBadPixels
, hithresh => 0.0045
, flickertimesteps => 1
, thresholdlabel => 'rate'
, withsearchbadcolumn => 'Y'
, hicolthresh => 0.00105
, backgroundrate => 0.00011
, loenergythresh => 0.14
, hienergythresh => 12.0
, narrowerthanpsf => 1.5
)
or return exception();
@extraopts = (
getnewbadpix => 'Y'
, badpixset => $newBadPixels
);
}
# Flag bad pixels
#
doCommand(
'badpix'
, eventset => $evFile1
, getuplnkbadpix => 'Y'
, getotherbadpix => 'Y'
, windowfilter => 'N'
, @extraopts
)
or return exception();
# Run epreject - reduce low-energy PN noise
# Imaging mode only
if ( $datamode eq "IMAGE.EL" ) {
doCommand(
'epreject'
, eventset => $evFile1
, set => $evFile0
, sigma => 4.0
# , noiseparameter => "0.98 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0"
, withoffsetmap => 'Y'
)
or return exception();
# run epnoise to reject soft X-ray noise in EPIC-pn
doCommand(
'epnoise'
, set => $evFile1
, identifynoisyframes => 'Y'
, noisecut => 2
, sigmacut => 3.0
, savemasks => 'N'
)
or return exception();
}
# Run epevents
#
my $evFile2 = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, ccd => $ccd
, content => 'Processed events'
, level => 2
);
info("Starting epevents....");
doCommand(
'epevents'
, eventset => $evFile1
, gainctiaccuracy => 2
, reemissionthresh => 0
, randomizeposition => 'Y'
, randomizeenergy => 'Y'
, withtempcorrection => 'Y'
, outset => $evFile2
)
or return exception();
info("Completed epevents....");
# Apply attitude correction
#
doCommand(
'attcalc'
, eventset => $evFile2
, refpointlabel => 'pnt'
, attitudelabel => 'ahf'
, withmedianpnt => 'Y'
, atthkset => $attFile
, -w => 10
)
or return exception();
# Filter events before final merging
my $evFile3 = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, ccd => $ccd
, content => 'Processed events'
, level => 3
);
# my $expression = '#XMMEA_EP && (PI>150 || PI<=150)';
my $expression = '(FLAG & 0x0fa0000)==0 && PI>150';
doCommand(
'evselect'
, table => "${evFile2}:EVENTS"
, withfilteredset => 'Y'
, filteredset => $evFile3
, destruct => 'Y'
, expression => $expression
, writedss => 'Y'
, updateexposure => 'Y'
, keepfilteroutput => 'Y'
)
or return exception();
# Append to list of processed event files
#
push( @finalEvents, $evFile3 );
# Accumulate expression for GTI filter
#
push( @gtiFilter, "(CCDNR==${ccd} && GTI($gtiFile1,TIME))" );
}
# Combine per-CCD event lists into single list
#
@extraopts = ();
my $tmpEvents;
if ( $datamode eq "IMAGE.EL" )
{
$tmpEvents = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Merged imaging events'
);
@extraopts = (
imagingset => $tmpEvents
, epnimgcolnames => [
qw(TIME RAWX RAWY DETX DETY X Y PHA PI FLAG PATTERN PAT_ID PAT_SEQ OFFSETX)
]
, epnimgcoltypes => [
qw(double int16 int16 int16 int16 int32 int32 int16 int16 int32 int8 int16 int8 int16)
]
);
}
elsif ( $datamode eq "TIME.EL" || $datamode eq "BURST.EL" )
{
$tmpEvents = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Merged timing events'
);
@extraopts = (
timingset => $tmpEvents
, epntimcolnames => [
qw(TIME RAWX RAWY DETX DETY X Y PHA PI FLAG PATTERN PAT_ID PAT_SEQ OFFSETX)
]
, epntimcoltypes => [
qw(double int16 int16 int16 int16 int32 int32 int16 int16 int32 int8 int16 int8 int16)
]
);
}
doCommand(
'evlistcomb'
, eventsets => [@finalEvents]
, instrument => 'epn'
, maintable => [qw(EVENTS OFFSETS)]
, othertables => [qw(BADPIX EXPOSURE DLIMAP)]
, @extraopts
)
or return exception();
# + + + Add epexposure call w/time randomisation, as requested by Carlos
doCommand(
'epexposure'
, eventsets => $tmpEvents
, screenexposure => 'Y'
, spatialexposure => 'N'
, randomizetime => 'Y'
) or return exception();
# Remove unwanted columns from EXPOSUnn extensions
my @badColumns = (qw()); # disabled, as epexposure does this already?
#my @badColumns = (qw(FRAMELIM NDSCLIN1 NDSCLIN2 NDSCLIN3 NDSCLIN4));
my @badObjects;
for ( my $ic = 1 ; $ic <= 12 ; $ic++ )
{
my $ext = sprintf 'EXPOSU%02d', $ic;
info("Checking for extension $ext");
if ( hasFITSExtension( file => $tmpEvents, extension => $ext ) )
{
my @columns = @badColumns;
push @badObjects
, map { $_ = "$tmpEvents:$ext:$_"; } @columns;
}
}
if (@badObjects)
{
doCommand( 'dsrm', objects => [@badObjects] )
or return exception();
}
# Fold in the HK gti files
#
my $hkGtiFileAligned = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Aligned HK GTIs'
);
if ( fileExists( file => $hkGTIFile ) )
{
doCommand(
"gtialign"
, gtitable => $hkGTIFile . ":STDGTI"
, eventset => $tmpEvents
, outset => $hkGtiFileAligned
)
or return exception();
foreach my $gtiFile (
findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'CCD GTI'
, level => 1
)
)
{
my $fileInf
= fileInfo( class => 'intermediate', name => $gtiFile );
my $gtiExt = sprintf 'STDGTI%02d', $$fileInf{'ccd'};
doCommand(
"gtimerge"
, tables => [ "$gtiFile", "$hkGtiFileAligned:$gtiExt" ]
, mergemode => 'and'
, withgtitable => 'N'
)
or return exception();
}
}
# Filter event list on GTI
#
my $gtiExpr = join( " || ", @gtiFilter );
my $finEvents;
if ( $datamode eq "IMAGE.EL" )
{
$finEvents = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'EPIC PN imaging mode event list'
);
}
elsif ( $datamode eq "TIME.EL" || $datamode eq "BURST.EL" )
{
$finEvents = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'EPIC timing mode event list'
);
}
doCommand(
'evselect'
, table => $tmpEvents
, filteredset => $finEvents
, withfilteredset => 'Y'
, expression => $gtiExpr
, keepfilteroutput => 'Y'
, writedss => 'Y'
, updateexposure => 'Y'
, destruct => 'true'
)
or return exception();
# Copy CIF into event list
#
my $cifFile = findFile(
class => 'product'
, content => 'Calibration index file'
);
doCommand(
'dscopyblock'
, blocks => "${cifFile}:CALINDEX"
, to => $finEvents
)
or return exception();
# Make flare GTI
#
if ( $datamode eq "IMAGE.EL" )
{
$timebinsize *= 2 if $mode eq "PrimeLargeWindow";
$timebinsize *= 35 if $mode eq "PrimeSmallWindow";
my $flareExpr = "($gtiExpr)";
if (@maskFilter)
{
$flareExpr .= " && (" . join( " || ", @maskFilter ) . ")";
}
# Save copy of bright pixel/CCD GTI filter for SFBLC processing
my $maskExpr = $flareExpr;
# Put in energy range for flare rejection
$flareExpr .= " && (PI IN [7000:15000])";
my $backRates = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Flare background timeseries'
);
my $flareGti = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'flare GTI'
);
doCommand(
'evselect'
, table => $finEvents
, maketimecolumn => 'Y'
, makeratecolumn => 'Y'
, timebinsize => $timebinsize
, withrateset => 'Y'
, rateset => $backRates
, writedss => 'Y'
, updateexposure => 'N'
, expression => $flareExpr
)
or return exception();
#
# 1 pixel ~ 17 arcsec**2 ~ 0.00474 arcmin**2, convert to per [ks * arcmin**2]
#
my $scalefactor = 0.0047355e-3; # unit = [arcmin**2/pixel ks/s]
# If rate file is empty, no rate extension is created. Trap this.
info("Number of good background pixels found: $goodPixels");
if ( hasFITSExtension( file => $backRates, extension => "RATE" )
&& $goodPixels > 100 )
{
# doCommand(
# 'tabcalc'
# , tables => "$backRates:RATE"
# , columntype => 'real64'
# , column => 'RATE'
# , expression =>
# "(COUNTS/($timebinsize*$goodPixels*$scalefactor))"
# )
# or return exception();
# Now subtract TSTART from the TIME values and leave the result in a new column T_ELAPSED:
my $start_time=readFITSKeyword(
file => $backRates
, extension => 'RATE'
, keyword => "TSTART"
);
doCommand(
'tabcalc'
, tables => "$backRates:RATE"
, columntype => 'real64'
, column => 'T_ELAPSED'
, expression => "TIME-$start_time" # Use the perl variable rather than the string '#TSTART' in case we decide later to use some other time as offset.
)
or return exception();
my $fplotCommandFile = newFile(class => 'intermediate'
, content => 'Flare bkg ts pco file'
, format => 'ASCII'
, extension => 'pco'
);
writeASCIIFile(
name => $fplotCommandFile
, text => [join("\n", "LA X TIME-$start_time (s)", 'LA G2 RATE (cts/s)', 'log y', 'PLOT', 'QUIT', '')]
);
my $flareBkgPs = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Flare bkg ts'
, format => 'PS'
);
doCommand(
'fplot'
, infile => $backRates
, xparm => 'T_ELAPSED'
, yparm => 'RATE'
, rows => '-'
, device => "${flareBkgPs}/cps"
, pltcmd => '@'."$fplotCommandFile"
)
or return exception();
#
# 1 pixel ~ 17 arcsec**2 ~ 0.00474 arcmin**2, convert to per [ks * arcmin**2]
#
# my $scalefactor = 0.0047355e-3; # unit = [arcmin**2/pixel ks/s]
#
# If rate file is empty, no rate extension is created. Trap this.
# info("Number for good background pixels found: $goodPixels");
# if ( hasFITSExtension( file => $backRates, extension => "RATE" )
# && $goodPixels > 100 )
# {
# doCommand(
# 'tabcalc'
# , tables => "$backRates:RATE"
# , columntype => 'real32'
# , column => 'RATE'
# , expression =>
# "(COUNTS/($timebinsize*$goodPixels*$scalefactor))"
# )
# or return exception();
# apply relative corrections to flare background timeseries intermediate
my $tmpBackRates = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Copy flare background timeseries'
);
doCommand(
'epiclccorr'
, srctslist => $backRates
, eventlist => $finEvents
, outset => $tmpBackRates
, applyabsolutecorrections => 'N'
, allcamera => 'Y'
) or return exception();
copyFile(
source => $tmpBackRates
, destination => $backRates
);
# Add bkgoptrate test call, progress towards flare optimised background subtraction (Mantis 73)
my $retval = 0;
my $optCmd= "bkgoptrate fracexpstyle=auto tssettabname=$backRates:RATE withendtime=N withlowercutoffcount=N withmintimeratio=N ignorelowcnttail=Y withstarttime=N -w 10";
info("RUN: $optCmd");
my $optRate = `$optCmd`;
chomp $optRate;
$retval = $?;
if ($retval ne 0) { return exception("Error $retval returned by bkgoptrate: $!"); }
info("OPTRATE: $optRate");
doCommand('tabgtigen', table => $backRates
, gtiset => $flareGti
, expression => "RATE<$optRate"
, mingtisize => $mingtisize
) or return exception();
# Insert source-free bkg flare light curve determination
# in-band source-free lightcurve generation for PN
# to include in pnEvents/sasflare8.5
info("==== Start in-band source-free flare lightcurve generation");
# Generate image for SFBLC analysis
my $SFBLCInputImage = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'SFBLC image'
, band => 8
, format => 'FITS'
);
# Find attitude GTI file
my $attGTIFile = findFile(
class => 'intermediate'
, instrument => 'all'
, content => 'attitude GTI'
);
# define inband image filter expression
my $inband_image_expr = "(RAWY > 12) && GTI($flareGti, TIME) && GTI($attGTIFile, TIME) && (( PATTERN <= 4 && PI in (500:1000] && (FLAG & 0x2fa002c) == 0) || (PATTERN<=4 && PI in (1000:7500] && (FLAG & 0x2fa0024) == 0))";
# make image (make use of existing flare GTI)
doCommand(
'evselect'
, table => $finEvents
, imageset => $SFBLCInputImage
, xcolumn => 'X'
, ycolumn => 'Y'
, imagebinning => 'binSize'
, ximagebinsize => 80
, yimagebinsize => 80
, expression => $inband_image_expr
, imagedatatype => 'Int32'
, withimagedatatype => 'true'
, writedss => 'true'
, updateexposure => 'true'
, keepfilteroutput => 'false'
) or return exception();
# check at least 1000 counts in input image before continuing to process
my $source_free_expr; # selection expression for source-free background light curve
# my $pn_timebinsize; downstream references to $pn_timebinsize changed to $timebinsize
my $sourceFreeBackRates;
my $SFBLCRegionFileFaint;
my $SFBLCRegionFileBright;
if (! &ModuleUtil::isImageLimit($SFBLCInputImage,'',1000.0)) {
# make vignetted exposure map
my $SFBLCExpmap = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'SFBLC exposure map'
, band => 8
, format => 'FITS'
);
doCommand(
'eexpmap'
, imageset => $SFBLCInputImage
, eventset => $finEvents
, attitudeset => $attFile
, attrebin => 2
, expimageset => $SFBLCExpmap
, pimin => 500
, pimax => 7500
, withvignetting => 'yes'
) or return exception();
# make non-vignetted exposure map
my $SFBLCNonVigExpmap = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'SFBLC non-vig exposure map'
, band => 8
, format => 'FITS'
);
doCommand(
'eexpmap'
, imageset => $SFBLCInputImage
, eventset => $finEvents
, attitudeset => $attFile
, attrebin => 2
, expimageset => $SFBLCNonVigExpmap
, pimin => 500
, pimax => 7500
, withvignetting => 'no'
) or return exception();
# make detection mask
my $SFBLCDetmask = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'SFBLC detmask'
, band => 8
, format => 'FITS'
);
doCommand(
'emask'
, detmaskset => $SFBLCDetmask
, expimageset => $SFBLCNonVigExpmap
, threshold1 => 0.5
, threshold2 => 1
);
# now run eboxdetect in local mode
# not trying to merge images
info("Running eboxdetect in local mode on EPIC PN, exposure $exp_id...");
# declare box-local source list
my $SFBLCBoxList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'SFBLC box-local source list'
, band => 8
, format => 'FITS'
);
# declare box-local source list (faint subset)
my $SFBLCBoxListFaint = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'SFBLC box-local source list faint subset'
, band => 8
, format => 'FITS'
);
# declare box-local source list (bright subset)
my $SFBLCBoxListBright = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'SFBLC box-local source list bright subset'
, band => 8
, format => 'FITS'
);
# generate box-local source list
doCommand(
'eboxdetect'
, imagesets => $SFBLCInputImage
, expimagesets => $SFBLCNonVigExpmap
, detmasksets => $SFBLCDetmask
, boxlistset => $SFBLCBoxList
, usemap => 'false'
, likemin => 5
, boxsize => 5
, withdetmask => 'true'
, nruns => 1
, ecf => "1.0"
, pimin => 500
, pimax => 7500
) or return exception();
# Now need to isolate ID_INT=ID_BAND=0 lines to a new file
# and convert to a region file
doCommand(
'fselect'
, infile => $SFBLCBoxList
, outfile => $SFBLCBoxListFaint
, expr => "(RATE <= 0.35)"
) or return exception();
doCommand(
'fselect'
, infile => $SFBLCBoxList
, outfile => $SFBLCBoxListBright
, expr => "(RATE > 0.35)"
) or return exception();
my $SFBLCTempFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'SFBLC temp file'
, band => 8
, format => 'FITS'
);
$SFBLCRegionFileFaint = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'SFBLC reg faint'
, band => 8
, format => 'FITS'
);
$SFBLCRegionFileBright = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'SFBLC reg bright'
, band => 8
, format => 'FITS'
);
my $lthresh = 50; # ML threshold for removing bright sources
my $rad = 60; # Fixed extraction rafius in arcsec (faint)
my $rad2 = 100; # Fixed extraction rafius in arcsec (bright)
my $region_select = "(ID_BAND == 1 && ID_INST == 1 && LIKE >= $lthresh && BOX_CTS > 0)";
doCommand(
'region'
, eventset => "${finEvents}:EVENTS"
, tempset => $SFBLCTempFile
, srclisttab => $SFBLCBoxListFaint
, expression => $region_select
, operationstyle => 'global'
, fixedradius => $rad
, bkgregionset => "${SFBLCRegionFileFaint}:REGION"
) or return exception();
doCommand(
'region'
, eventset => "${finEvents}:EVENTS"
, tempset => $SFBLCTempFile
, srclisttab => $SFBLCBoxListBright
, expression => $region_select
, operationstyle => 'global'
, fixedradius => $rad2
, bkgregionset => "${SFBLCRegionFileBright}:REGION"
) or return exception();
# Generate temporary evaluation image
my $SFBLCEvalImage = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'SFBLC temp evaluation image'
, band => 8
, format => 'FITS'
);
# Include regions in eval selection expression
my $eval_image_expr = $inband_image_expr . " && (region($SFBLCRegionFileFaint,X,Y) && region($SFBLCRegionFileBright,X,Y))";
doCommand(
'evselect'
, table => $finEvents
, imageset => $SFBLCEvalImage
, xcolumn => 'X'
, ycolumn => 'Y'
, withimageset => 'Y'
, imagebinning => 'binSize'
, ximagebinsize => 80
, yimagebinsize => 80
, expression => $eval_image_expr
, imagedatatype => 'Int32'
, withimagedatatype => 'true'
, writedss => 'true'
, updateexposure => 'true'
, keepfilteroutput => 'false'
) or return exception();
# Now, finally should be able to use the region file on the event list to
# generate the light curve for flare filtering
# PN
# need to establish the exposure ID and other parts of the name
# needs the CCD GTI and bright pixel masks
# Am following the PN light curve creation steps which look more
# complicated/extensive than for MOS
# info("Creating CCD/GTI filters for EPIC PN");
# my @gtiexpr;
# foreach my $gtiFile (
# findFile(
# class => 'intermediate'
# , instrument => thisInstrument
# , exp_id => $exp_id
# , content => 'CCD GTI'
# , level => 1
# )
# )
# {
# my $fileInf
# = fileInfo( class => 'intermediate', name => $gtiFile );
# my $ccd = $$fileInf{'ccd'};
# push @gtiexpr, "( CCDNR==$ccd && GTI($gtiFile, TIME))";
# };
# my $gtiexpr_string = "(" . join(" || ", @gtiexpr) . ")";
#
# info("Creating Bright Pixel Masks for EPIC PN");
#
# my @brtexpr;
#
# foreach my $brtmask (
# findFile(
# class => 'intermediate'
# , instrument => thisInstrument
# , exp_id => $exp_id
# , content => 'Bright pixels mask'
# )
# )
# {
#
# my $fileInf
# = fileInfo( class => 'intermediate', name => $brtmask );
# my $ccd = $$fileInf{'ccd'};
# push @brtexpr, "( CCDNR == $ccd && MASK($brtmask,1,1,RAWX,RAWY))";
# };
# my $brtexpr_string = "(" . join(" || ", @brtexpr) . ")";
# merge faint and bright region files prior to generating source free light curve
my $SFBLCRegionFileMerged = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'SFBLC reg merged'
, band => 8
, format => 'FITS'
);
doCommand(
'fmerge'
, infiles => [$SFBLCRegionFileBright, $SFBLCRegionFileFaint]
, outfile => $SFBLCRegionFileMerged
, columns => "SHAPE,X,Y,R,ROTANG,COMPONENT"
, mextname => 'REGION'
, history => 'Y'
, copyprime => 'Y'
) or return exception();
doCommand(
'fparkey'
, value => 'pos'
, fitsfile => "${SFBLCRegionFileMerged}[REGION]"
, keyword => "MTYPE1"
, add => 'Y'
) or return exception();
doCommand(
'fparkey'
, value => 'X,Y'
, fitsfile => "${SFBLCRegionFileMerged}[REGION]"
, keyword => "MFORM1"
, add => 'Y'
) or return exception();
# make in-band source-free light curve
info("Creating source-free in band light curve for PN");
# $source_free_expr = $maskExpr." && (PATTERN <=4) && (PI IN [500:7500]) && ((FLAG & 0xfffffef) == 0) && (region($SFBLCRegionFileFaint,X,Y) && region($SFBLCRegionFileBright,X,Y))";
$source_free_expr = $maskExpr." && (PATTERN <=4) && (PI IN [500:7500]) && ((FLAG & 0xfffffef) == 0) && (region($SFBLCRegionFileMerged,X,Y))";
} else {
# extraction expression for sparse images (omit region select)
info("Creating in band light curve for PN (sparse image < 1000 counts)");
$source_free_expr = $maskExpr." && (PATTERN <=4) && (PI IN [500:7500]) && ((FLAG & 0xfffffef) == 0)";
}
# Generate source-free background timeseries as product
# $pn_timebinsize = 10;
$sourceFreeBackRates = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'EPIC flare background timeseries'
);
doCommand(
'evselect'
, expression => $source_free_expr
, maketimecolumn => 'Y'
, makeratecolumn => 'Y'
, rateset => $sourceFreeBackRates
, table => $finEvents
, timebinsize => $timebinsize
, updateexposure => 'N'
, withrateset => 'Y'
, writedss => 'Y'
) or return exception();
# add light curve rate and T_ELAPSED columns. The RATE column definition
# differs from current pipeline value because there the cut is made by
# rate per unit sky area. Here we don't care about unit area so we can
# just divide the counts by the bin time
# info("Adding column for RATE, using PN time binsize $pn_timebinsize");
#
# doCommand(
# 'tabcalc'
# , column => 'RATE'
# , columntype => 'real64'
# , expression => "COUNTS/${pn_timebinsize}"
# , tables => "${sourceFreeBackRates}:RATE"
# ) or return exception();
# Observation start time should already be determined in pnEvents
# Apply relative corrections to output product timeseries
my $tmpSourceFreeBackRates = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Copy EPIC flare background timeseries'
);
doCommand(
'epiclccorr'
, srctslist => $sourceFreeBackRates
, eventlist => $finEvents
, outset => $tmpSourceFreeBackRates
, applyabsolutecorrections => 'N'
, allcamera => 'Y'
) or return exception();
copyFile(
source => $tmpSourceFreeBackRates
, destination => $sourceFreeBackRates
);
# Generate T_ELAPSED column
info("Adding column for T_ELAPSED, using observation start time $start_time");
doCommand(
'tabcalc'
, column => 'T_ELAPSED'
, columntype => 'real64'
, expression => "TIME-$start_time"
, tables => "${sourceFreeBackRates}:RATE"
) or return exception();
# Run bkgoptrate - optimum cut threshold
my $snDiagFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Flare GTI SN diagnostic'
, band => 8
, format => 'FITS'
);
my $retval = 0;
my $rcut_pn = `bkgoptrate -V 0 -w 0 tssettabname=${sourceFreeBackRates}:RATE dooutputsntab=yes ignorelowcnttail=yes snsettabname=${snDiagFile}:SN_TO_BKGCUT`;
info("CMD: bkgoptrate -V 0 -w 0 tssettabname=${sourceFreeBackRates}:RATE dooutputsntab=yes ignorelowcnttail=yes snsettabname=${snDiagFile}:SN_TO_BKGCUT");
$retval = $?;
if ($retval ne 0) { return exception("Error $retval returned by bkgoptrate: $!"); }
my $optFlareGti;
if ($rcut_pn eq '') {
info("bkgoptrate call returns no value, skipping optimised flare GTI generation...");
} else {
if (fileExists( file => $sourceFreeBackRates ) ) {
info("bkgoptrate returns optimum cut threshold: $rcut_pn");
# append SN_TO_BKGCUT table to FBKTSR product
doCommand(
'ftappend'
, infile => "${snDiagFile}[SN_TO_BKGCUT]"
, outfile => $sourceFreeBackRates
) or return exception();
# add FLCUTTHR header value
doCommand(
'addattribute'
, set => "${sourceFreeBackRates}:RATE"
, attributename => 'FLCUTTHR'
, attributetype => 'real'
, attributecomment => '"\"Optimised flare cut threshold\""'
, realvalue => $rcut_pn
) or return exception();
# cut threshold limited to 100
# if ($rcut_pn > 100) {
# $rcut_pn = 100;
# info("Optimum cut threshold limited to 100");
# }
# Now generate the optimised flare GTI files
$optFlareGti = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => "Optimised flare GTI"
, band => 8
, format => 'FITS'
);
info("Creating GTI files with this threshold");
doCommand(
'tabgtigen'
, table => $sourceFreeBackRates
, gtiset => $optFlareGti
, expression => "RATE < ${rcut_pn}"
, mingtisize => 100
) or return exception();
}
}
# make plot of background light curve
info("Make full in-band, source-free light curve plots");
my $sourceFreeLightCurvePS = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => "In-band source-free light curve"
, band => 8
, format => 'PS'
);
my $sourceFreeLCCommandFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'In-band sf lc pco file'
, format => 'ASCII'
, extension => 'pco'
);
# Cutoff threshold label string
my $srcut = sprintf("%.3f",$rcut_pn);
# Write PLT command file
writeASCIIFile(
name => $sourceFreeLCCommandFile
, text => [join("\n","LA X TIME-$start_time (s)",'LA G2 RATE (cts/s)','log y','lwid 5',"LA 1 P 0 $srcut LI 0 1 CO 4 \"$srcut\"",'PLOT','QUIT','')]
);
doCommand(
'fplot'
, infile => $sourceFreeBackRates
, xparm => 'T_ELAPSED'
, yparm => 'RATE'
, rows => '-'
, device => "${sourceFreeLightCurvePS}/CPS"
, pltcmd => '@'."$sourceFreeLCCommandFile"
) or return exception();
# Get extremal values of calculated S/N
my $snBg = readFITSColumn(
file => $snDiagFile
, extension => 'SN_TO_BKGCUT'
, column => 'SN_RATIO'
) or return exception();
my $snMin = min(@$snBg);
my $snMax = max(@$snBg);
# Generate plot of S/N
my $SNRCurvePS = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'In-band sf Signal to noise plot'
, band => 8
, format => 'PS'
);
my $SNRCommandFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'In-band sf Signal to noise pco file'
, format => 'ASCII'
, extension => 'pco'
);
writeASCIIFile(
name => $SNRCommandFile
,text => [join("\n","LA X Background rate cut (cts/s)",'LA G2 Signal-to-noise','lwid 5',"LA 1 Pos $srcut $snMax LI -90 1.0 CO 4 \"$srcut\"",'PLOT','QUIT','')]
);
doCommand(
'fplot'
, infile => $snDiagFile
, xparm => 'BKGRATECUT'
, yparm => 'SN_RATIO'
, rows => '-'
, device => "${SNRCurvePS}/CPS"
, pltcmd => '@'."$SNRCommandFile"
) or return exception();
my $tempPS = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'In-band sf temp plot file'
, band => 8
, format => 'PS'
);
my $combiPS = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'In-band sf combi plot file'
, band => 8
, format => 'PS'
);
$retval=0;
my $cli_1 = `gs -sDEVICE=pswrite -sOutputFile=${tempPS} -dNOPAUSE -dBATCH $sourceFreeLightCurvePS $SNRCurvePS`;
info("CMD: gs -sDEVICE=pswrite -sOutputFile=${tempPS} -dNOPAUSE -dBATCH $sourceFreeLightCurvePS $SNRCurvePS");
$retval = $?;
if ($retval ne 0) { return exception("Error returned by gs: $!"); }
$retval=0;
my $cli_2 = `psnup -s0.6 -2up ${tempPS} > ${combiPS}`;
info("CMD: psnup -s0.6 -2up ${tempPS} > ${combiPS}");
$retval=$?;
if ($retval ne 0) { return exception("Error returned by psnup: $!"); }
# Generate optimised flare background PDF
my $flareBkgPdf = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'EPIC flare background timeseries'
, format => 'PDF'
);
PStoPDF(
source => $combiPS
, destination => $flareBkgPdf
)
or return exception();
# Regenerate light curve but filter by new GTIs and plot to see if they disappear
info("TESTING: Recreating light curves with only valid GTIs and regions left in");
my $optBackRates = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Optimised flare bkg time series'
, band => 8
, format => 'FITS'
);
my $test_ts_expr = $maskExpr." && (PATTERN <=4) && GTI($attGTIFile, TIME) && (PI in [500:7500]) && ((FLAG & 0xfffffef) == 0)";
if (fileExists( file => $SFBLCRegionFileFaint )) {
$test_ts_expr .= " && region($SFBLCRegionFileFaint,X,Y)";
}
if (fileExists( file => $SFBLCRegionFileBright )) {
$test_ts_expr .= " && region($SFBLCRegionFileBright,X,Y)";
}
if (fileExists( file => $optFlareGti )) {
$test_ts_expr .= " && GTI($optFlareGti, TIME)";
}
doCommand(
'evselect'
, expression => $test_ts_expr
, maketimecolumn => "Y"
, makeratecolumn => 'Y'
, rateset => $optBackRates
, table => $finEvents
, timebinsize => $timebinsize
, updateexposure => 'N'
, withrateset => 'Y'
, writedss => 'Y'
) or return exception();
# doCommand(
# 'tabcalc'
# , column => 'RATE'
# , columntype => 'real64'
# , expression => "COUNTS/${pn_timebinsize}"
# , tables => "${optBackRates}:RATE"
# ) or return exception();
# Apply relative corrections to test GTI timeseries
my $tmpOptBackRates = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Copy optimised flare bkg time series'
, band => 8
, format => 'FITS'
);
# Require valid nonzero flare GTI, and non-zero RATE column for final epiclccorr call
if( (numberFITSRows(file=>$optFlareGti, extension=>'STDGTI') > 0) &&
!ModuleUtil::isColumnEmpty($optBackRates,"RATE") )
{
doCommand(
'epiclccorr'
, srctslist => $optBackRates
, eventlist => $finEvents
, outset => $tmpOptBackRates
, applyabsolutecorrections => 'N'
, allcamera => 'Y'
) or return exception();
copyFile(
source => $tmpOptBackRates
, destination => $optBackRates
);
doCommand(
'tabcalc'
, column => 'T_ELAPSED'
, columntype => 'real64'
, expression => "TIME-$start_time"
, tables => "${optBackRates}:RATE"
) or return exception();
my $OptSourceFreeLightCurvePS = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Opt in-band source-free light curve'
, band => 8
, format => 'PS'
);
doCommand(
'fplot'
, infile => $optBackRates
, xparm => 'T_ELAPSED'
, yparm => 'RATE'
, rows => '-'
, device => "${OptSourceFreeLightCurvePS}/CPS"
, pltcmd => '@'."$sourceFreeLCCommandFile"
) or return exception();
}
}
}
# Everything OK
return success();
}
1;