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.46';
$version = $VERSION;
$author='Dean Hinshaw,Richard West,Duncan John Fyfe,Ian Stewart,Duncan Law-Green,Eduardo Ojero,Ed Chapin,Jose V Perea';
$date='2014-06-25';
#
# ChangeLog
# =========
#
# Version 3.46 - 2014-06-25 (JVP)
# ------------
#
# + Plot for 'Flare GTI SN diagnostic' produced from extension [SN_TO_BKGCUT] of product 'EPIC flare background timeseries'
# Filename labels in plots are final products now, instead intermediate files.
#
# Version 3.45 - 2013-08-26 (EC)
# ------------
#
# + RATE column from evselect must be non-zero for final epiclccorr call
#
# Version 3.44 - 2013-06-11 (EC)
# ------------
#
# + Merge changes from SOC version.
#
# version 3.431 - 2012-11-22 EO
# ----------------------------
#
# - Added execution of command attcalc after execution of evselect (line 1225)
# if condition (fileExists(file=>$tmpimgEvents) is fulfilled (line 1217),
# to compute the full image of a mosaic observation.
#
# version 3.43 - 2012-11-05 DLG
# -----------------------------
#
# - merge faint and bright region files for generation of source-free background light curve
#
# version 3.42 - 2012-11-05 DLG
# -----------------------------
#
# - require valid flare GTI for final epiclccorr call (crashes on empty GTI)
#
# version 3.41 - 2012-10-25 DLG
# -----------------------------
#
# - Modify plotting of SN output graph
#
# version 3.40 - 2012-10-24 DLG
# -----------------------------
#
# - evselect makeratecolumn=Y, disable subsequent tabcalc
#
# version 3.39 - 2012-10-22 DLG
# -----------------------------
#
# - all allcamera = Y to epiclccorr call
#
# version 3.38 - 2012-10-17 DLG
# -----------------------------
#
# - check for existence of SFBLC region files in evselect
#
# version 3.37 - 2012-10-09 DLG
# -----------------------------
#
# - set writedss=Y on remaining SSP evselect calls
# - remove writedss param from epiclccorr
#
# version 3.36 - 2012-10-08 DLG
# -----------------------------
#
# - try writedss=Y on epiclccorr call
# + move tablcalc T_ELAPSED after epiclccorr
#
# version 3.35 - 2012-10-01 DLG
# -----------------------------
#
# - reenable epiclccorr, specify outset as temp file
#
# version 3.34 - 2012-10-01 DLG
# -----------------------------
#
# - disable epiclccorr call pending GTI crash fix
#
# version 3.33 - 2012-09-27 DLG
# -----------------------------
#
# + correct epiclccorr param to 'applyabsolutecorrections'
#
# version 3.32 - 2012-09-26 DLG
# -----------------------------
#
# + Apply relative time series corrections using epiclccorr()
#
# version 3.31 - 2012-09-26 DLG
# -----------------------------
#
# + Undo change 3.30 (moved to EPICSourceProducts)
# + correct evselect logic for region files
# + modify RATE column to real64
#
# version 3.30 - 2012-09-04 DLG
# -----------------------------
#
# + copy region file extensions to FBKTSR product
#
# version 3.29 - 2012-09-04 DLG
# -----------------------------
#
# + Alter time binning on flare background timeseries to 26s (from 52s)
#
# version 3.28 - 2012-09-04 DLG
# -----------------------------
#
# + Add source regions to evselect for Optimised flare bkg time series
#
# version 3.27 - 2012-08-22 DLG
# -----------------------------
#
# + Generate SFBLC and SNR dual plot, with optimised SN marked
#
# version 3.26 - 2012-08-15 DLG
# -----------------------------
#
# + append SN_TO_BKGCUT table on background time series product
#
# version 3.25 - 2012-07-27 DLG
# -----------------------------
#
# + modify SFBLC region code: source exclusion radius to be dependent on source counts
#
# version 3.24 - 2012-07-23 DLG
# -----------------------------
#
# + disable 100cts maximum on $rcut_mos
#
# version 3.23 - 2012-07-10 DLG
# -----------------------------
#
# + reinstate bkgoptrate call for bad pixel filtering
#
# version 3.22 - 2012-07-09 DLG
# -----------------------------
#
# + add ignorelowcnttail=yes parameter for bkgoptrate calls
#
# Version 3.21 - 2011-07-02 DLG
# -----------------------------
#
# + Refactor code for sparse image handling
#
# Version 3.20 - 2011-07-02 DLG
# -----------------------------
#
# + Trap null return from bkgoptrate, skip optimised GTI generation
#
# Version 3.19 - 2011-06-25 DLG
# -----------------------------
#
# + Remove global attitude GTI filter constraint fron in-band SF light curve
#
# Version 3.18 - 2011-06-25 DLG
# -----------------------------
#
# + Add FLCUTTHR header keyword to opt flare light curve
#
# Version 3.17 - 2011-06-22 DLG
# -----------------------------
#
# + Do not perform source-free background analysis for sparse images < 1000 counts
#
# Version 3.16 - 2011-06-21 DLG
# -----------------------------
#
# + Produce optimised flare background time series FITS and PDF as product (FBKTSR)
#
# Version 3.15 - 2011-06-20 DLG
# -----------------------------
#
# + Clamp maximum threshold to 100 for flare GTI generation
#
# Version 3.14 - 2011-06-10 DLG
# -----------------------------
#
# + Disable bkgoptrate call for bad pixel threshold
#
# Version 3.13 - 2011-06-07 DLG
# -----------------------------
#
# + Add SN diagnostic for source-free bkgoptrate calc
#
# Version 3.12 - 2011-05-25 DLG
# -----------------------------
#
# + Add source-free background flare lightcurve generation
#
# Version 3.11 - 2011-04-06 DLG
# -----------------------------
#
# + Add second bkgoptrate test, 1st for bright pixel check, 2nd for optimised GTI
#
# Version 3.10 - 2011-03-07 DLG
# ------------
#
# + Add bkgoptrate test, optimised flare background filtering (Mantis 73)
#
# Version 3.091 - 2013-01-09 EC
# -------------
#
# + Add calls to manualintervention to allow skipping of CCDs
#
# Version 3.09 - 2009-10-05 DLG
# ------------
#
# + Plot EPIC background flare light curves in log space (Mantis 52)
#
# Version 3.08 - 2009-09-20 DLG
# ------------
#
# + Added counting cycle report check, set $mingtisize accordingly
#
# Version 3.07 - 2009-07-28 DLG
# ------------
#
# + Set mingtisize = 0 for testing
#
# Version 3.06a - 2009-06-19 DLG
# -------------
#
# + Set ignore() if SAS_ATTITUDE environment variable is 'RAF' for slew processing
#
# Version 3.06 - 2005-11-16 IMS
# ------------
#
# + Incorporated DJF's shortened intermediate file names.
#
# Version 3.05 - 2005-10-31 IMS
# ------------
#
# + Parameter --getotherbadpix of badpix call changed from 'no' to 'yes' on advice from JB (should get rid of bright pixels in ccd 4, mos1.)
#
# Version 3.04 - 2005-10-14 IMS
# ------------
#
# + TSTART obtained from flare time series header and subtracted from the times of the flare PDF plot.
# + Flare plot Y axis now RATE not COUNTS. Calculation of RATE is thus moved forward a bit.
# + Added 'PStoPDF', 'writeASCIIFile' to the list of methods called from ModuleResources.
#
# Version 3.03 - 2005-09-30 RGW
# ------------
#
# + OFFSETS table is propagated by evlistcomb
#
# version 3.02 - 2005-08-12 RGW
# ------------
#
# + trimmed filenames in flare plotting to work around FPLOT filename length
# restriction
#
# version 3.01 - 2005-08-09 IMS
# ------------
#
# + PDF product now made from flare bkg time series.
# + Calls to fmedian removed because embadpixfind>2.0 now calculates this internally. Parameter --medianset of embadpixfind has likewise been removed.
#
# version 3.00 - 2005-05-09 DJF
# ------------
#
# + Add explicit perl module headers. Previously these were supplied
# at compile time. This will make debugging and extending the modules
# through additional perl libraries easier.
#
# Version 2.33 - 2004-03-19 (DJF)
# ------------
#
# + Added 'FLICKERING' and 'ON_BADOFFSET' to flareFlags and eventlistFlags.
#
# Version 2.32 - 2004-03-17 (DJF)
# ------------
#
# + Added Explicit use of withimageset=>'N' for evselect.
#
# Version 2.31 - 2004-02-12 (DJF)
# ------------
#
# + Use new PCMS function gtiTimeThreshold() to return threshold
#
# Version 2.30 - 2004-01-14 (DJF)
# ------------
#
# + Use tabgitgen parameter mingtisize rather than evselect to remove
# short intervals from bad pixel and flare background GTI.
#
# Version 2.29 - 2003-12-09 (DJF)
# ------------
#
# + Adapted doCommand to use anonymous lists for list parameters
#
# Version 2.28 - 2003-07-01 (DJF)
# ------------
#
# + #XMMEA_EM selection had been applied to wrong evselect. This caused
# removal of GATTI events before flare background construction.
# Try again. #XMMEA_EM selection applied to final event list
#
# Version 2.27 - 2003-07-01 (DJF)
# ------------
#
# + XMMEA_EM selection still not working, more testing.
#
# Version 2.26 - 2003-07-01 (DJF)
# ------------
#
# + Fixed use of #XMMEA_EM flag selection in good event filtering.
# From (FLAG & #XMMEA_EM)==0 to (#XMMEA_EM)==0
#
# Version 2.25 - 2003-07-01 (DJF)
# ------------
#
# + Corrected good event filtering. Had missed out the ==0 from the flag filter.
#
# Version 2.24 - 2003-06-30 (DJF)
# ------------
#
# + Added emosdatamodes => 'IMAGING' to first evlistcomb call.
#
# Version 2.23 - 2003-04-11 (DJF)
# ------------
#
# + Fixed command line to evlistcomb (instrument= 'emos' rather than thisInstrument)
#
# Version 2.22 - 2003-04-08 (DJF)
# ------------
#
# + Added Jean Ballet's newer recipy for improved flare screening
#
# Version 2.21 - 2002-07-01 (DJF)
# ------------
#
# + Changed image creation flags to XMMEA_EM
#
# Version 2.20 - 2002-03-21 (DH)
# ------------
#
# + Bug fix for version 2.19.
#
# Version 2.19 - 2002-03-20 (DH)
# ------------
#
# + Put in removal of short (<100s) flare GTIs for bad pixel
# finding.
# + Use offset variance files in first call to emevents.
# + Add bitmask fitlering to final event list, so it makes it in
# the data subspace.
#
# Version 2.18 - 2002-03-13 (DH)
# ------------
#
# + Fix bugs in identifying and handling of CCDs in low gain mode.
#
# Version 2.17 - 2002-03-11 (DH)
# ------------
#
# + Bug fix: was using wrong file name for output flare gti file.
#
# Version 2.16 - 2002-03-05 (DH)
# ------------
#
# + Fix error in background rate calculation. Do rate calculation
# based on the ccd area available, rather than the number of CCDs.
#
# Version 2.15 - 2002-02-20 (DH)
# ------------
#
# + Back out change in 2.14, as the MOS software is not yet ready to
# handle this.
#
# Version 2.14 - 2002-02-19 (DH)
# ------------
#
# + Add OFFSETS to the list of tables for evlistcomb to propagate.
#
# Version 2.13 - 2002-02-18 (DH)
# ------------
#
# + Explicitly list all parameters for the fmedian ftool call.
# + More bug fixes
#
# Version 2.12 - 2002-02-15 (DH)
# ------------
#
# + Bug fixes for 2.11 .
#
# Version 2.11 - 2002-02-07 (DH)
# ------------
#
# + New method for specifying bitmask filters for event lists, in order
# to make it more clear which flags are being used.
# + Add in the UNDERSHOOT flag, for event list and flare screening.
# + The flare background timeseries is now generated as a rates file, not
# a histogram. The RATE column is normalized to 1/[ks*acrmin**2], and
# the GTI cut is made at 2.0 .
# + Check for GAIN_CCD keyword in frame file. If set to 'LOW', do not
# process events for this ccd.
#
# Version 2.10 - 2002-01-08 (DH)
# ------------
#
# + Test for XMMEA_22 in event list before making flare histograms.
#
# Version 2.09 - 2001-11-14 (DH)
# ------------
#
# + Fix bug in flare GTI application introduced in 2.07.
#
# Version 2.08 - 2001-11-06 (DH)
# ------------
#
# + Slight change to flare GTI test to more exactly mimic MakeImage.pm .
#
# Version 2.07 - 2001-10-31 (DH)
# ------------
#
# + Test flare GTI as in MakeImage.pm to see if it is above a time threshold before using it.
#
# Version 2.06 - 2001-08-20 (DJF)
# ------------
#
# + Added test for empty GTI when filtering out flares before searching for bright pixels.
#
# Version 2.05 - 2001-08-15 (DJF)
# ------------
#
# + Added test for column 'HISTO' to first call to tabgtigen.
#
# Version 2.04 - 2001-08-14 (DJF)
# ------------
#
# + Remove flagtruncatede1=N flag from first call to emevents (to fix XMMEA_22 failure fo evselect)
# + Add flag findbright=N to first call to embadpixfind as suggesteg by Jean Ballet.
#
# Version 2.03 - 2001-08-13 (DJF)
# ------------
#
# + Another fix for histogram intermadiate to product change.
#
# Version 2.02 - 2001-08-13 (DJF)
# ------------
#
# + Fix to flare histogram intermadiate to product change.
#
#
# Version 2.01 - 2001-08-09 (DJF)
# ------------
#
# + Brought up to date with changes to previous versions of MOSEvents.
# + Changed Flare Histogram from intermediate to final product.
# - Removed supurfluous findFile at end of module v1.13
#
#############################################
## Changes Applied to v1.10 before #####
## implimentation of v 2.0.1 #####
##
## Version 1.13 - 2001-06-07 (DH)
## ------------
##
## + Take out changes in version 1.12 and move test for
## empty event list to MakeImage module.
##
## Version 1.12 - 2001-06-07 (DH)
## ------------
##
## + Set ignore flag if final imaging event list is empty.
##
## Version 1.11 - 2001-06-07 (DH)
## ------------
##
## + Do not make flare gti if flare histogram is empty.
##
############################################
#
# Version 2.00 - 2001-05-29 (DH)
# ------------
#
# + Complete re-work of bad pixel finding using the new embadpixfind package.
#
# Version 1.10 - 2001-05-21 (DH)
# ------------
#
# + Changes made as instructed by Jean Ballet:
# - changes to rejection of e3 and e4 events due to task changes
# - Add call to emenergy before badpixfind in Imaging mode only
# - Change to emenergy options for TIME and CTIME modes
#
# Version 1.09 - 2001-03-16 (DH)
# ------------
#
# + Print out version number in performAction() for
# tracking purposes.
#
# Version 1.08 - 2001-01-11 (DH)
# ------------
#
# + Remove GlobalHK from ignore rule.
#
# + Add InstrumentHK to start rule.
#
# Version 1.07 - 2000-11-23 (DH)
# ------------
#
# + First production version.
#
# Declare list of instruments this module is interested in
@instList=qw(emos1 emos2);
# Number of streams
sub numberOfStreams {
return numberOfExposures();
}
sub evaluateRules {
# If upstream module has been ignored, so if this one
return ignore() if ignored(module => 'ExpChooser',
instrument => thisInstrument,
stream => thisStream);
# ignore for slew processing
return ignore() if ($ENV{'PCMS_ISSLEW'});
# Start conditions
start()
if complete(module => 'ExpChooser',
instrument => thisInstrument,
stream => thisStream)
and
complete(module => 'InstrumentHK',
instrument => thisInstrument,
stream => 1)
and
complete(module => 'GlobalHK',
instrument => 'all',
stream => 1);
}
my $flareFlags = ['UNDERSHOOT'
, 'BAD_E3E4'
, 'ON_BADROW'
, 'OUTSIDE_THRESHOLDS'
, 'OUT_OF_CCD_WINDOW'
, 'ON_BADPIX'
, 'COSMIC_RAY'
, 'IN_BAD_FRAME'
, 'OUT_OF_FOV'
, 'FLICKERING'
, 'ON_BADOFFSET'
];
my $eventlistFlags =['UNDERSHOOT'
, 'BAD_E3E4'
, 'ON_BADROW'
, 'OUTSIDE_THRESHOLDS'
, 'OUT_OF_CCD_WINDOW'
, 'ON_BADPIX'
, 'COSMIC_RAY'
, 'IN_BAD_FRAME'
, 'FLICKERING'
, 'ON_BADOFFSET'
];
sub performAction {
info("Module version number: $version");
# Misc. variables
my @extraopts;
# Work out which exposure this stream maps on to
my $exp_id=exposureID(instrument => thisInstrument
,stream => thisStream
);
info("Exposure ID: $exp_id");
#
my $filter=getExposureProperty(exp_id => $exp_id
,name => 'filter'
);
my $mode=getExposureProperty(exp_id => $exp_id
,name => 'mode'
);
info("Filter: $filter");
info("Instrument mode: $mode");
# Find event files for this exposure
my @evFiles=findFile(class => 'ODF'
,instrument => thisInstrument
,exp_id => $exp_id
,content => '*events'
);
# Return immediately if there are no event files
if ($#evFiles<0)
{
ppd("000005");
return success();
}
my @evFiles2;
# Find auxiliary file for this exposure in the ODF
my $auxFile=findFile(class => 'ODF'
,instrument => thisInstrument
,exp_id => $exp_id
,content => 'Auxiliary data'
,required => 'true'
);
# Find attitude HK file
my $attFile=findFile(class => 'product'
,instrument => 'all'
,content => 'Attitude time series'
,required => 'true'
);
# Find HK GTI file
my $hkGTIFile=findFile(class => 'intermediate'
,instrument => thisInstrument
,content => 'HK GTI'
);
# Set mingtisize to 100 by default
my $mingtisize = 100;
# Check for counting cycle report
# reset mingtisize if required
my @ccxfiles = findFile(class => 'ODF'
,instrument => thisInstrument
,exp_id => $exp_id
,content => 'Counting cycle report'
);
unless (@ccxfiles) {
info("No counting cycle report found. Setting mingtisize=0");
$mingtisize = 0;
}
# Loop through each event file
my @filter;
my $fovArea = 0;
foreach my $evFile0 (@evFiles)
{
# Find CCD and node numbers
my $ccd=readFITSKeyword(file => $evFile0,
extension => 2,
keyword => "CCDID"
);
# manualintervention will check for keyword/keyval
# pairs provided to the "--ignore" command-line option
# of "sequence" that match the supplied "instrument"
# and "exp_id" values, respectively. Although
# manualintervention was originally provided to skip
# entire exposures, there is no reason we can't use it
# here, provided that the keyword for skipping CCDs of
# an instrument is not identical to a pure instrument
# name which would cause things to be skipped earlier
# in ExpChooser. To skip CCDs we append "ccd" to the
# instrument name.
# e.g.,
# sequence --ignore="emos1ccd=4,7"
# will skip CCDs 4 and 7 of MOS1.
if (manualintervention( instrument => thisInstrument . "ccd", exposure => $ccd ))
{
info("Manual Intervention. Ignoring CCD $ccd");
next;
}
my $node=readFITSKeyword(file => $evFile0,
extension => 2,
keyword => "CCDNODE"
);
my $datamode=readFITSKeyword(file => $evFile0,
extension => 2,
keyword => "DATATYPE"
);
info("First pass of events for CCD $ccd, node $node, mode=$datamode");
# Run emframes
#
my $frameFile=newFile(class => 'intermediate',
instrument => thisInstrument,
exp_id => $exp_id,
node => $node, 'ccd' => $ccd,
content => 'Rectified frames'
);
my $gtiFile1=newFile(class => 'intermediate',
instrument => thisInstrument,
exp_id => $exp_id,
node => $node,
ccd => $ccd,
content => 'CCD GTI',
level => 1
);
my $evFile1=newFile(class => 'intermediate',
instrument => thisInstrument,
exp_id => $exp_id,
node => $node,
ccd => $ccd,
content => 'Processed events',
level => 1
);
@extraopts=();
@extraopts=(flagbadtimes => 'yes', ingtiset => $hkGTIFile)
if $hkGTIFile
;
doCommand('emframes',
auxiliaryset => $auxFile,
frameset => $frameFile,
withodfeventset => 'Y',
odfeventset => $evFile0,
neweventset => 'Y',
outeventset => $evFile1,
writegtiset => 'Y',
outgtiset => $gtiFile1,
@extraopts
) or return exception();
if( hasFITSKeyword(file=>$frameFile, extension=>'FRAMES', keyword=>'GAIN_CCD')
&& readFITSKeyword(file=>$frameFile, extension=>'FRAMES', keyword=>'GAIN_CCD') =~ /^LOW/
){
info("CCD $ccd is in low gain mode. Will not process events for this ccd.");
ppd(id => "000007" , ccd => $ccd , node => $node );
next;
}
push @evFiles2, $evFile0;
#
@extraopts=();
if( $datamode eq "IMAGE.EL" || $datamode eq "RIMAGE.EL" )
{
my $evFile2=newFile(class => 'intermediate',
instrument => thisInstrument,
exp_id => $exp_id,
node => $node,
ccd => $ccd,
content => 'Processed events',
level => 2
);
@extraopts=();
@extraopts=( rejectbade3 => 'N' )
unless ($datamode eq "IMAGE.EL");
# Check if an offset-variance map is available for this
# instrument/CCD/node, if so add to the arguments for 'emevents'
#
my @ovfile=findFile(class => 'ODF',
instrument => thisInstrument,
ccd => $ccd,
node => $node,
content => 'Offset/variance mode events'
);
push(@extraopts, withoffvarsets => 'Y', offvarsets => join(" ",@ovfile))
if @ovfile;
ppd(id => "000006" , ccd => $ccd , node => $node) if @ovfile;
doCommand('emevents',
odfeventset => $evFile1,
eventset => $evFile2,
analysepatterns => 'N',
flagbadpixels => 'N',
splitdiagonals => 'N',
withframeset => 'Y',
frameset => $frameFile,
randomizeposition => 'N',
@extraopts)
or return exception();
doCommand('emenergy',
ineventset => $evFile2,
correctcti => 'N',
correctgain => 'N',
randomizeenergy => 'N')
or return exception();
$fovArea += readFITSKeyword(file => $evFile2
,extension => 'EVENTS'
,keyword => 'IN_FOV'
) if ( hasFITSKeyword(file => $evFile2
,extension => 'EVENTS'
,keyword => 'XMMEA_22'
)
);
}
}
# Make flare gti for use with bad pixel searching
my $hasBadpixGti = 0;
my $bpFlareGti;
my $backgrndBinsize = 26;
my $GTI_time_threshold=gtiTimeThreshold();
my $fovAreaThreshold = 10;
if( $fovArea > $fovAreaThreshold )
{
ppd(id => "000008" , data => { fovarea => $fovArea , threshold => $fovAreaThreshold } );
my $bpEvFile=newFile(class => 'intermediate',
instrument => thisInstrument,
exp_id => $exp_id,
content => 'Bad pix find event list'
);
my @outerEvFiles=findFile(class => 'intermediate'
,instrument => thisInstrument
,exp_id => $exp_id
,content => 'Processed events'
,level => 2
);
doCommand('evlistcomb'
,eventsets => [@outerEvFiles]
,imagingset => $bpEvFile
,instrument => 'emos'
,emosdatamodes => 'IMAGING'
);
$bpFlareGti=newFile(class => 'intermediate',
instrument => thisInstrument,
exp_id => $exp_id,
content => 'Bad pix flare GTI');
my $histoFile=newFile(class => 'intermediate',
instrument => thisInstrument,
exp_id => $exp_id,
content => 'Flare rates file',
level => 1);
if( hasFITSKeyword(file => $bpEvFile, extension => 'EVENTS', keyword => 'XMMEA_22') )
{
my $flareBitmask = getEventBitmask(file=>$bpEvFile, masks=>$flareFlags);
ppd(id => "000009" , data => { flarebitmask => $flareBitmask } );
doCommand('evselect'
,table => $bpEvFile
,writedss => 'N'
,updateexposure => 'N'
,withrateset => 'Y'
,rateset => $histoFile
,maketimecolumn => 'Y'
,makeratecolumn => 'Y'
,timebinsize => $backgrndBinsize
,timecolumn => 'TIME'
,expression => "(PATTERN==0) && #XMMEA_22 && ((FLAG & $flareBitmask) == 0)"
,withimageset => 'N'
) or return exception();
if( numberFITSRows(file=>$histoFile, extension=>'RATE') > 0 )
{
# 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();
# Added to make the full image of a mosaic observation.
doCommand('attcalc',
,eventset => $imgEvents,
,refpointlabel => 'pnt',
,withmedianpnt => 'true',
,attitudelabel => 'ahf',
,atthkset => $attFile,
,calctlmax => 'yes')
or return exception();
# 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();
doCommand(
'fplot'
, infile => "${sourceFreeBackRates}[SN_TO_BKGCUT]"
, xparm => 'BKGRATECUT'
, yparm => 'SN_RATIO'
, rows => '-'
, device => "${SNRCurvePS}/CPS"
, pltcmd => '@'."$SNRCommandFile"
) or return exception();
my $tempPS = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'In-band sf temp plot file'
, band => 8
, format => 'PS'
);
my $combiPS = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'In-band sf combi plot file'
, band => 8
, format => 'PS'
);
$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;