package ExpDetect;
use strict;
use English;
use Carp;
use vars
qw(@ISA $VERSION $name $author $date $version @instList $numberOfStreams);
@ISA = qw(Module);
use ModuleResources;
use ModuleUtil;
# Declare identity, version, author, date, etc.
$name = __PACKAGE__;
$VERSION = '2.421';
$version = $VERSION;
$author = 'Richard West,Dean Hinshaw,Duncan John Fyfe,Richard Saxton,Ian Stewart,Duncan Law-Green,Ed Chapin';
$date = '2013-10-21';
#
# ChangeLog
# =========
# version 2.421 - 2013-10-21 EC
# ------------
#
# + An additional tweak regarding ep_updatevig: If emldetect doesn't
# find any sources, ep_updatevig will fail with an exception. Trap
# this case by checking for a non-empty source list first.
#
# version 2.42 - 2013-08-14 EC
# ------------
#
# + Temporarily add call to ep_updatevig (local package provided by
# Iris Traulsen) to fix incorrect vignetting factor columns output by
# emldetect. Once emldetect is replaced we can revert to the old
# behaviour
#
# version 2.41 - 2012-11-18 DLG
# ------------
#
# + (temporary?) disable removal of LE CCDs from $emTagFilteredList
# emtaglenoise does not appear to detect LE noise consistently
#
# version 2.40 - 2012-08-29 DLG
# ------------
#
# + eootepileupmask ratethreshold param is instrument-dependent
#
# version 2.39 - 2012-08-20 DLG
# ------------
#
# + modify dpssflag call to version 11, pass ootergamasklist, minstlist params
#
# version 2.38 - 2012-07-03 DLG
# ------------
#
# + correct dpssflag to --withflag10=no in case of no valid OOTE masks
#
# version 2.37 - 2012-06-27 DLG
# ------------
#
# + in case of no per-camera OOTE masks, dpssflag runs without otrgamaskset parameter set
#
# version 2.36 - 2012-05-02 DLG
# ------------
#
# + Add threshold=10 to XID emldetect call
#
# version 2.35 - 2012-05-01 DLG
# ------------
#
# + Add psfmodel = ellbeta to both emldetect calls
#
# version 2.34 - 2012-04-30 DLG
# ------------
#
# + Add forrgaspike=yes to eootepileupmask call
#
# version 2.33 - 2012-04-25 DLG
# ------------
#
# + Amalgamate call to dpssflag from SrcMerge to prevent flag overwrite
#
# version 2.32 - 2012-04-23 DLG
# ------------
#
# + Check for existence of $filteredList before starting first half of MOS LE noise filtering
#
# version 2.31 - 2012-04-19 DLG
# ------------
#
# + Add dpssflag to mark OOTE sources from pileup mask
#
# version 2.30 - 2012-04-16 DLG
# ------------
#
# + Start implementation of OOT masking (eootpileupmask)
#
# version 2.29 - 2012-03-23 DLG
# ------------
#
# + Explicitly write merged exposure maps as EXPMAP products
#
# version 2.28 - 2012-03-21 DLG
# ------------
#
# + separate emtaglenoise and evselect calls, due to bug which wipes DSS header keywords
#
# version 2.27 - 2012-03-13 DLG
# ------------
#
# + Dataflow correction by LH. @filtImList tracks MOS LE filtered images
#
# version 2.26 - 2012-03-12 DLG
# ------------
#
# + Correction: stage 3 eexpmap call should use LE Filtered image
#
# version 2.26 - 2012-03-07 DLG
# ------------
#
# + Input image merged event list is band-8 to support emtaglenoise
#
# version 2.25 - 2012-02-16 DLG
# ------------
#
# + Comtinue MOS LE noise correction
#
# version 2.24 - 2011-12-06 DLG
# ------------
#
# + Start MOS LE noise correction
#
# version 2.23 - 2010-06-10 DLG
# ------------
#
# + PN dark columns fix (Mantis #xxx)
#
# version 2.22 - 2009-12-10 DLG
# ------------
#
# + add psfmodel = 'ellbeta' for PSF model testing in psfdev
#
# version 2.21 - 2009-12-01 DLG
# ------------
#
# + emdetect changes for PSF testing and development
#
# version 2.20 - 2006-08-16 DJF
# ------------
#
# + ONTIME and EXPOSURE fix moved to ExpDetectPrep where it belongs.
#
# version 2.19 - 2006-08-15 DJF
# ------------
#
# + Fix ONTIME and EXPOSURE in cross-camera merged images.
#
# version 2.18 - 2006-07-13 DJF
# ------------
#
# + Saving image information for later emosaic use was using epic instrument
# rather than individual instrument name in filename
#
# version 2.17 - 2006-07-11 DJF
# ------------
#
# + Detection mask promoted to a full product. Without the detection mask there is no
# obvious means of determining why no sources were detected with a non-zero exposure.
#
# version 2.16 - 2006-03-22 DJF
# ------------
#
# + Corrected parameters to keepSafe function
#
# version 2.15 - 2006-03-08 DJF
# ------------
#
# + Changed eboxdetect likemin to 5.
# + esensemap likemin set to that of X as advised by Y
#
# version 2.14 - 2006-02-21 DJF
# ------------
#
# + Exposure maps which are supplied to the band 8 imweightadd image are
# written to a file to be used with emosaic in ImageMerge module.
# Using the imweightadd output looses the indidual image keywords.
#
# version 2.13 - 2005-11-29 IMS
# ------------
# + Value of emask parameter --threshold1 changed from 0.05 to 0.5
#
# version 2.12 - 2005-11-04 DJF
# ------------
# + Changed 'merged bkg' to 'merged background' for band 8 images to make it consistent with ImageMerge module.
#
# version 2.11 - 2005-11-03 IMS
# ------------
#
# + Fixed a bug in the sub isInList.
# + The second call to emldetect now proceeds only if there are the same number of image files as ecfs.
#
# version 2.10 - 2005-11-02 IMS
# ------------
#
# + Non-vig exposure map now has a shorter content string (thus shorter name).
# + No longer make PNG products for bkg maps in bands 1-5.
# + Band 8 fits bkg map changed from product to intermediate.
# + Band 8 bkg map PNG product abolished.
#
# version 2.09 - 2005-10-31 IMS
# ------------
#
# + Updated task parameter values according to scripts sent by GL on 23 Oct 05.
# + XID (merged-exposure) bkg maps now being made.
# + Final emldetect call now uses XID images, exposure maps etc made in the ExpDetectPrep module.
# + 'Observation' changed to 'merged' in the content string for bkg maps.
# + Sensitivity maps now being made.
#
# version 2.08 - 2005-10-20 IMS
# ------------
#
# + Keywords listing the constituent exposure ids are now written to box and mldetect source list headers.
# + Same now also written to background maps.
# + Added GIFtoPNG to list of used ModuleResources subroutines.
# + The non-vignetted rather than the vignetted exp map is now sent to emask.
#
# version 2.07 - 2005-10-19 RGW
# ------------
#
# + shortened name of intermediate EPIC observation background map files (80-character filename limits nibbles the privates yet again).
#
# version 2.06 - 2005-09-23 RGW
# ------------
#
# + deal gracefully with cases where merged maps don't exist for an instrument
#
# version 2.05 - 2005-09-21 IMS
# ------------
#
# + Some debugging done, so it will run.
#
# version 2.04 - 2005-09-20 IMS
# ------------
#
# + FITS background map products are now also being made for band 8.
# + PNG background maps now made for bands 1-5, 8.
# + Slight relocation of the generation of the band 1-5 background map names.
#
# version 2.03 - 2005-09-14 IMS
# ------------
#
# + The first emldetect output (contains info for bands 1-5) is now made the product, rather than the second (which just contains band 9 info).
#
# version 2.02 - 2005-08-16 RGW
# ------------
#
# + handling of lists of mask images corrected
# + changed spline map to a product rather than intermediate
#
# version 2.01 - 2005-07-12 IMS
# ------------
#
# + All instrument-specific parts made non-instrument-specific, since we are now doing a single EPIC source detection.
# + $numberOfStreams=1 instead of =numberOfExposures().
# + evaluateRules changed to reflect new prior module ExpDetectPrep.
# + $exp_id and $mode are no longer necessary (this processing is now done in ExpDetectPrep).
# + Images, exp maps and filtered source lists which have been merged in ExpDetectPrep are now read as inputs rather than single-exposure ones. (If there's only 1 exposure per inst+filter combination, ExpDetectPrep copies the single-exposure image/expmap/list to the 'merged' name.)
# + Non-empty $filteredList is now necessary for push of $rawim, $expim etc into lists.
# + Filter is now read from the band-0 merged source list for each instrument.
# + Det mask made for each instrument and added to list @maskImList.
# + Parameter settings are those recommended by GL on 15 Jun 2005.
# + XID now added in a separate, final call to emldetect.
# + Non-vignetted, merged exposure map is now also used (supplied to esplinemap).
#
# + ***>>> NOTE! <<<*** esensmap invocation temporarily commented out. What to do? Alternatives: (i) put in a loop over insts; (ii) single call to esensitivity to make single epic sens map.
#
# 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
# + Standardized date format (ccyy-mm-dd)
# + Standardized author list to use comma separator
# + Make use of perl $VERSION magic. Now $Version = version = ''
#
# Version 1.38 - 2004-03-03
# ------------
#
# + Moved energy conversion factors into ModuleResources.
# + Moved hardness bands to ModuleResources.
# + Changed energyBands() match changes in ModuleResources.
#
# Version 1.37 - 2004-02-19
# ------------
#
# + Fixed evaluateRules cross-instrument dependance
#
# Version 1.36 - 2004-01-14
# ------------
#
# + Changed dependance on MakeImage to MakeMOSImage and MakePNImage
#
# Version 1.35 - 2003-12-09
# ------------
#
# + Adapted doCommand to use anonymous lists for list parameters
#
# Version 1.34 - 2003-06-13 (RDS)
# ------------
#
# + New ECFs for latest QE files
#
# Version 1.33 - 2002-09-16 (DJF)
# ------------
#
# + Remove code for testing spline node variations
#
# Version 1.32 - 2002-08-27 (DJF)
# ------------
#
# + Fix background image list for testing of esplinemap
#
# Version 1.31 - 2002-08-23 (DJF)
# ------------
#
# + Try different node values in esplinemap
#
# Version 1.30 - 2002-08-16 (DJF)
# ------------
#
# + Change to esplinemap parameters. Use scut = 0.001 rather than default
# + Added call to esplinemap and eboxdetect for different number of nodes.
#
# Version 1.29 - 2002-06-10 (DJF)
# ------------
#
# + Incorporating catalogue pipeline improvements into the production pipeline
#
# (See patch versions below for further details)
#
# Version 1.28.12 - 2002-05-08 (DH)
# ---------------
#
# + Remove eposcorr call for ml exposure source lists.
#
# Version 1.28.11 - 2002-05-02 (DH)
# ---------------
#
# + Add pimin and pimax parameters to the esplinemap call.
#
# Version 1.28.10 - 2002-04-30 (DH)
# ---------------
#
# + Fix typo.
#
# Version 1.28.9 - 2002-04-26 (DH)
# --------------
#
# + Make sure implot only tries to plot FITS files.
#
# Version 1.28.8 - 2002-04-26 (DH)
# --------------
#
# + Put oot options back in esplinemap.
#
# Version 1.28.7 - 2002-04-26 (DH)
# --------------
#
# + For now remove oot options from esplinempap, unitl bug
# is fixed.
#
# Version 1.28.6 - 2002-04-25 (DH)
# --------------
#
# + Bug fix.
#
# Version 1.28.5 - 2002-04-24 (DH)
# --------------
#
# + Shorten name of intermediate image event lists to
# accomodate a bug in esplinemap.
#
# Version 1.28.4 - 2002-04-19 (DH)
# --------------
#
# + Add ml list and exp map graphics products.
# + Add new feature in esplinemap for handling oot events.
#
# Version 1.28.3 - 2002-03-29 (DH)
# --------------
#
# + Correct spelling of eposcorr.
#
# Version 1.28.2 - 2002-03-29 (DH)
# --------------
#
# + Bug fixes for 1.28.2.
# + For extended source detection, need to do a call to eboxdetect
# as well as emldetect.
# + Add in call to eposscorr to correct source positions
#
# Version 1.28.1 - 2002-03-28 (DH)
# --------------
#
# + Call to esky2det and extended source detection on xid band.
#
# Version 1.28 - 2002-03-28 (DH)
# ------------
#
# + BKGDSCRN keyword now always written (either T or F). Comment
# fixed.
#
# Version 1.27 - 2002-02-25 (DH)
# ------------
#
# + Put pn band 1 image back into source detection.
#
# Version 1.26 - 2002-02-08 (DH)
# ------------
#
# + Fix bug in esplinemap call. The parameter idband should be the band
# reference number, rather than the actual band number itself.
#
# Version 1.25 - 2002-02-08 (DH)
# ------------
#
# + Bug fixes for 1.24.
#
# Version 1.24 - 2002-02-08 (DH)
# ------------
#
# + Add keyword BKGDSCRN to products if flare GTI was applied to images.
#
# Version 1.23 - 2002-01-07 (DH)
# ------------
#
# + Don't run emldetect if the box map source list is empty.
#
# Version 1.22 - 2001-10-30 (DH)
# ------------
#
# + Put in new ECF parameters, as per SSC-LUX-TN-0059, issue 2.0. Includes
# values for the Open filter position.
#
# Version 1.21 - 2001-10-09 (DH)
# ------------
#
# + Fix mistake in changing parameter name hrdef in eboxdetect.
#
# Version 1.20 - 2001-10-08 (DH)
# ------------
#
# + Replace hrdef and xiddef with new parmeters which specify
# the instrument (hrm1def, hrm2def, ...);
# + Change value of likemin parameter from 10 to 8. Change mlmin
# parameter to 10, whereas is was 16 or 18 depending on the number
# of bands.
#
# Version 1.19 - 2001-10-04 (DH)
# ------------
#
# + Revert back to old ECFs of version 1.17 .
#
# Version 1.18 - 2001-10-01 (DH)
# ------------
#
# + Correct MOS think filter ECFs. Closes SSC-SPR-2573.
#
# Version 1.17 - 2001-08-30 (DFJ)
# ------------
#
# + corrected ecf instrument name, pn -> epn
#
# Version 1.16 - 2001-08-21 (DFJ)
# ------------
#
# + Correct some miss typed energy conversion factors
#
# Version 1.15 - 2001-08-15 (DFJ)
# ------------
#
# + Apply correction to use of new xidecf parameter.
#
# Version 1.14 - 2001-08-14 (DFJ)
# ------------
#
# + Corrected energy conversion factors (see SSC-LUX-TM-0059).
# + Added XID band conversion factors (withecf) to emldetect
#
# Version 1.13 - 2001-06-04 (DH)
# ------------
#
# + Change nsplinenode of esplinemap from 10 to 16.
#
# Version 1.12 - 2001-05-30 (DH)
# ------------
#
# + Change mlmin value for emldetect task from 10 to 16 for 4 bands
# and 18 for 5 bands.
#
# Version 1.11 - 2001-05-21 (DH)
# ------------
#
# + Change esplinemap parameter withexpimage to false, and remove
# expimageset paramter.
#
# Version 1.10 - 2001-05-02 (DH)
# ------------
#
# + Add values for hrdef and efc parameters to eboxdetect. Use
# same values as for emldetect.
# + Remove pointresponse from emldetect, which has been removed
# from the latest version.
# + Add pimin and pimax parameters to emldetect and eboxdetect.
# Use the new energyBands function to obtain the needed values.
# + Add new xiddef parameter to emldetect, which gives which input
# images correspond to the xidbands.
# + Fix SSC-SPR-2361. The code was not using the proper xid bands
# for making the sensitivity maps.
#
# Version 1.09 - 2001-03-20 (DH)
# ------------
#
# + New version of esplinemap, which no longer has an idinst
# parameter. So take it out.
#
# Version 1.08 - 2001-03-19 (DH)
# ------------
#
# + Do not perform source detection for PN small window modes.
#
# Version 1.07 - 2001-03-16 (DH)
# ------------
#
# + Fix bug which was introduced in 1.05/1.06 which
# caused hrdef and idinst parameters not to set
# in esplinemap.
# + Print out version number in performAction() for
# tracking purposes.
#
# Version 1.06 - 2001-03-08 (DH)
# ------------
#
# + Set withxidband=true in emldetect.
# + Put energy conversion factors into emldetect.
# + Add in new 'hrdef' parameter to emldetect.
#
# Version 1.05 - 2001-03-06 (DH)
# ------------
#
# + Fix bug introduced in 1.04 which caused wrong energy
# bands to be used to make pn sensitivity maps.
#
# + Change calling sequence for new version of esplinemap:
# Set the idband and idinst parameters.
# Take the default setting for scut, whose value has been
# redefined.
#
# Version 1.04 - 2001-01-29 (DH)
# ------------
#
# + Remove energy band 1 from list of images used for source
# detection in the pn.
#
# Version 1.01 - 2000-12-11 (DH)
# ------------
#
# + First production version.
# Declare list of instruments this module is interested in
@instList = qw(epic);
#my $energyBands = [qw(1 2 3 4 5)];
my $xidBands = {'epn' => [qw(2 3 4)], 'emos1' => [qw(2 3 4)], 'emos2' => [qw(2 3 4)] };
### This list of xid bands ought to be obtained from a subroutine call because it is also needed in ExpDetectPrep. IMS 2005-10-31.
# Number of streams
sub numberOfStreams
{
1;
}
#
sub evaluateRules
{
return ignore()
if (
allIgnored( module => 'ExpDetectPrep' )
);
start()
if (
allComplete( module => 'ExpDetectPrep' )
);
}
#
sub performAction
{
info("Module version number: $version");
#
my (
@rawImList, @expImList, @flatExpImList, @bkgImList
, @ebandlo, @ebandhi, @imEventList, @maskImList
, @instList, @maskImLongList
, @bandLongList, @hardnessLongList, @ecfList
, @xidEbandlo, @xidEbandhi, @xidEcfList
, %band8constituents, %band9constituents, %expidKwds
, $emTagFilteredList, @filtImList, @lebadccds
);
# Set which bands to use for source searching
my @bandList = ( 1, 2, 3, 4, 5 );
my @hardnessList = qw(1 2 2 3 3 4);
# define cameras for LE noise correction
my %runLENoiseFilter = (
'epn' => 'n',
'emos1' => 'y',
'emos2' => 'y'
);
# Find raw images, exposure maps for this exposure, in band order
foreach my $inst (qw(epn emos1 emos2))
{
my $filteredList = findFile(
class => 'intermediate'
, instrument => $inst
, band => 8
, content => 'image merged event list'
);
# Start OOT masking ------------------------------------------------------------
# check for image merged event list
if ($filteredList) {
my $raw8im = findFile(
class => 'intermediate'
, instrument => $inst
, band => 8
, content => 'merged image'
, format => 'FITS'
, required => 'true'
);
my $attFile = findFile(
class => 'product'
, instrument => 'all'
, content => 'Attitude time series'
, required => 'true'
);
my $attRebin = newFile(
class => 'intermediate'
, instrument => 'all'
, content => 'Rebinned attitude time series'
, format => 'FITS'
);
doCommand(
'attbin'
, eventset => $filteredList
, attsource => 'atthk'
, atthkset => $attFile
, binnedattset => $attRebin
, maxdelta => 1.0
);
my $pileupMask = newFile(
class => 'intermediate'
, instrument => $inst
, content => 'Pileup mask'
, band => 8
, format => 'FITS'
);
# eootepileupmask flagging threshold is instrument dependent
my $ratethreshold;
if ($inst eq 'epn') {
$ratethreshold = 0.008; # PN value
} else {
$ratethreshold = 0.004; # MOS value
};
doCommand(
'eootepileupmask'
, ratethreshold => $ratethreshold
, forrgaspike => 'yes'
, eventset => $filteredList
, maskset => $pileupMask
, outputstyle => 'sky'
, templateset => $raw8im
, binnedattset => $attRebin
);
}
# Start LE noise filtering -----------------------------------------------------
if ($runLENoiseFilter{$inst} eq 'y' && -e $filteredList) {
info("Band-8 $inst event list found. Starting $inst LE noise filtering");
$emTagFilteredList = newFile(
class => 'intermediate'
, instrument => $inst
, band => 8
, content => 'LE noise tagged merged event list'
, format => 'FITS'
);
# Step 1: Flag bad CCDs --------------------------------------------------------
# emtaglenoise generates LENOISnn=1 keywords for CCDs affected by LE noise
doCommand(
'emtaglenoise'
, eventset => $filteredList
, filterbadccds => 'no'
) or return exception();
# Loop over CCDs and get value of LENOISnn keyword
@lebadccds = ();
foreach my $ccdnum ( 1 .. 12 ) {
my $lenoisekey = sprintf("LENOIS%02d", $ccdnum);
my $lenoisekeyval = readFITSKeyword(
file => $filteredList
, extension => 'PRIMARY'
, keyword => $lenoisekey
);
# info("DEBUG: Keyword $lenoisekey = $lenoisekeyval");
# DISABLE REMOVAL OF LE CCDS
# awaiting fix to emtaglenoise -- does not appear to detect LE noise consistently
#
# if ( $lenoisekeyval eq "1" ) {
# push (@lebadccds, "( CCDNR != $ccdnum )");
# }
}
# Construct LE noise selection expression
my $leexpr = join ( " && ", @lebadccds );
# LE filtered temp image
my $dummyimg = newFile(
class => 'intermediate'
, instrument => $inst
, band => 8
, content => 'LE filtered temp image'
, format => 'FITS'
);
# Run evselect to remove affected CCDs from band-8 merged event list
doCommand(
'evselect'
, expression => $leexpr
, filteredset => $emTagFilteredList
, imagebinning => "binSize"
, imagedatatype => "Int32"
, imageset => $dummyimg
, keepfilteroutput => "true"
, table => $filteredList
, updateexposure => "true"
, withimagedatatype => "true"
, writedss => "true"
, xcolumn => "X"
, ximagebinsize => 80
, ycolumn => "Y"
, yimagebinsize => 80
);
}
# Loop over bands -------------------------------------------------------------
my @xidBandsThisInst = @{$xidBands->{$inst}};
my ($filter, $maskIm);
foreach my $band (@bandList)
{
my $rawim = findFile(
class => 'intermediate'
, instrument => $inst
, band => $band
, content => 'merged image'
);
my $expim = findFile(
class => 'intermediate'
, instrument => $inst
, band => $band
, content => 'merged exposure map'
);
my $flatExpIm = findFile(
class => 'intermediate'
, instrument => $inst
, band => $band
, content => 'merged non-vig exp map'
);
if ( $rawim && $expim && $flatExpIm && $filteredList )
{
# Continue LE noise filtering -----------------------------------------------------
if ($runLENoiseFilter{$inst} eq 'y') {
info("Continuing $inst LE noise filtering");
# Step 2a: Create LE unfiltered image -------------------------------------------
# Get PI limits for this energy band
my $pilo = energyBand($band, 0);
my $pihi = energyBand($band, 1);
# declare LE filtered event list
my $LEFilteredList = newFile(
class => 'intermediate'
, instrument => $inst
, band => $band
, content => 'LE noise filtered image event list'
, format => 'FITS'
);
# declare LE unfiltered image
my $LEUnfilteredImage = newFile(
class => 'intermediate'
, instrument => $inst
, band => $band
, content => 'LE noise unfiltered image'
, format => 'FITS'
);
# Generate image
doCommand(
'evselect'
, table => $filteredList.":EVENTS"
, imageset => $LEUnfilteredImage
, xcolumn => "X"
, ycolumn => "Y"
, imagebinning => "binSize"
, ximagebinsize => 80
, yimagebinsize => 80
, squarepixels => 'true'
, expression => "(PI in $pilo:$pihi)"
, imagedatatype => "Int32"
, withimagedatatype => "true"
, writedss => "true"
, updateexposure => "true"
, keepfilteroutput => "y"
, filteredset => $LEFilteredList
) or return exception();
# Step 2b: Create LE filtered image -----------------------------------------
# declare LE filtered image
my $LEFilteredImage = newFile(
class => 'intermediate'
, instrument => $inst
, band => $band
, content => 'LE noise filtered image'
, format => 'FITS'
);
# Generate image
doCommand(
'evselect'
, table => $emTagFilteredList.":EVENTS"
, imageset => $LEFilteredImage
, xcolumn => "X"
, ycolumn => "Y"
, imagebinning => "binSize"
, ximagebinsize => 80
, yimagebinsize => 80
, squarepixels => 'true'
, expression => "(PI in $pilo:$pihi)"
, imagedatatype => "Int32"
, withimagedatatype => "true"
, writedss => "true"
, updateexposure => "true"
, keepfilteroutput => "y"
, filteredset => $LEFilteredList
) or return exception();
# Step 3: Run eexpmap to create exposure map ----------------------------------
# Need both with and without vignetting
# Find spacecraft attitude file
my $attFile = findFile(
class => 'product'
, instrument => 'all'
, content => 'Attitude time series'
, required => 'true'
);
# Define vignetted exposure map
my $LEVigExpmap = newFile(
class => 'intermediate'
, instrument => $inst
, band => $band
, content => 'LE noise filtered vignetted exposure map'
, format => 'FITS'
);
doCommand(
'eexpmap'
, imageset => $LEFilteredImage
, eventset => $emTagFilteredList
, attitudeset => $attFile
, attrebin => 2
, expimageset => $LEVigExpmap
, pimin => $pilo
, pimax => $pihi
, withvignetting => 'yes'
) or return exception();
my $LEUnVigExpmap = newFile(
class => 'intermediate'
, instrument => $inst
, band => $band
, content => 'LE noise filtered unvignetted exposure map'
, format => 'FITS'
);
doCommand(
'eexpmap'
, imageset => $LEFilteredImage
, eventset => $emTagFilteredList
, attitudeset => $attFile
, attrebin => 2
, expimageset => $LEUnVigExpmap
, pimin => $pilo
, pimax => $pihi
, withvignetting => 'no'
) or return exception();
# Generate product merged exposure map
my $expFile = newFile(
class => 'product'
, instrument => $inst
, band => $band
, content => 'EPIC merged exposure map'
, format => 'FITS'
);
copyFile(
source => $LEVigExpmap
, destination => $expFile
);
# Create detection mask for this exposure
info("Creating detection mask for $inst, band $band -- LE noise flagging enabled");
$maskIm = newFile(
class => 'product'
, instrument => $inst
, band => $band
, content => 'EPIC DETECTION MASK'
);
doCommand(
'emask'
, expimageset => $LEUnVigExpmap
, threshold1 => 0.5
, threshold2 => 1.0
, detmaskset => $maskIm
) or return exception();
push ( @maskImList, $maskIm );
push( @filtImList, $LEFilteredImage );
# Set LE noise filtered image as input for step 4
$rawim = $LEUnfilteredImage;
# Set LE vignetted exposure map as input for step 4
$expim = $LEVigExpmap;
# Set LE unvignetted exposure map as input for step 5
$flatExpIm = $LEUnVigExpmap;
} else {
# Generate product merged exposure map
my $expFile = newFile(
class => 'product'
, instrument => $inst
, band => $band
, content => 'EPIC merged exposure map'
, format => 'FITS'
);
copyFile(
source => $expim
, destination => $expFile
);
# Create detection mask for this exposure
info("Creating detection mask for $inst, band $band");
$maskIm = newFile(
class => 'product'
, instrument => $inst
, band => $band
, content => 'EPIC DETECTION MASK'
);
doCommand(
'emask'
, expimageset => $flatExpIm
, threshold1 => 0.5
, threshold2 => 1.0
, detmaskset => $maskIm
) or return exception();
push( @maskImList, $maskIm );
# Default to raw unfiltered image
push ( @filtImList, $rawim );
}
# =================================
if ( !$filter ) # ie, if this is the first band
{
$filter = readFITSKeyword(
file => $filteredList
, extension => 'PRIMARY'
, keyword => 'FILTER'
);
unless ( $filter =~ /(Thin|Medium|Thick|Open)/ )
{
info("Could not match filter for energy conversion factor");
return exception();
}
$filter = $1;
push(
@xidEcfList
, energyConversionFactor(
$inst, $filter, 9
)
);
push( @xidEbandlo, energyBand( 9, 0 ) );
push( @xidEbandhi, energyBand( 9, 1 ) );
# Read the constituent exposure IDs of the band 1 merged image:
my @expidKwdValues = ();
my $exited_normally = 0;
foreach my $i (1 .. 99) # why the arbitrary 99? Because I don't like infinite loops.
{
my $kwd_str = 'EXPID'.sprintf("%2.2d", $i);
my $kwd_value = readFITSKeyword(
'file' => $rawim
, 'extension' => 'PRIMARY'
, 'keyword' => $kwd_str
);
if ($kwd_value)
{
push( @expidKwdValues, $kwd_value);
}
else
{
$exited_normally = 1;
last;
}
}
if ($exited_normally)
{
$expidKwds{$inst} = [@expidKwdValues];
}
else
{
info("Bug when reading expid keywords.");
}
} # end if !$filter (ie, if $band==1)
my $bkgFile = newFile(
class => 'product'
, instrument => $inst
, band => $band
, content => 'EPIC merged background map'
);
push( @bkgImList, $bkgFile );
push( @{$band8constituents{$inst}}, $bkgFile);
push( @{$band9constituents{$inst}}, $bkgFile) if (&isInList($band, @xidBandsThisInst));
# info("DEBUG: Pushing $maskIm onto maskImLongList");
push( @maskImLongList, $maskIm );
push( @rawImList, $rawim );
push( @expImList, $expim );
push( @flatExpImList, $flatExpIm );
push(
@ecfList
, energyConversionFactor(
$inst, $filter, $band
)
);
push( @ebandlo, energyBand( $band, 0 ) );
push( @ebandhi, energyBand( $band, 1 ) );
push( @imEventList, $filteredList );
push( @instList, $inst );
push( @bandLongList, $band );
}
else
{
info(
"Cannot find all required images for instrument $inst, energy band $band. This combination will be ignored."
);
}
} # end loop over bands
} # end loop over insts
return success() unless @rawImList;
# Populate Band-1 detmask list
# iterate over instruments
my @band1MaskList;
foreach my $inst (qw(epn emos1 emos2))
{
my $band1Mask = findFile(
class => 'product'
, instrument => $inst
, band => 1
, content => 'EPIC DETECTION MASK'
);
push( @band1MaskList, $band1Mask );
}
my $flareScreened;
if (
hasFITSKeyword(
file => $rawImList[0]
, extension => 'PRIMARY'
, keyword => 'BKGDSCRN'
)
)
{
$flareScreened = readFITSKeyword(
file => $rawImList[0]
, extension => 'PRIMARY'
, keyword => 'BKGDSCRN'
);
undef $flareScreened if $flareScreened eq 'F';
}
# Step 4: Do first box-detection stage -------------------------------------
#
my $boxList = newFile(
class => 'product'
, instrument => thisInstrument
, content => 'EPIC OBSERVATION box-local source list'
);
doCommand(
'eboxdetect'
, imagesets => [@rawImList]
, expimagesets => [@expImList]
, detmasksets => [@band1MaskList]
, boxlistset => $boxList
, usemap => 'false'
, likemin => 5
, boxsize => 5
, withdetmask => 'true'
, withexpimage => 'true'
, nruns => 1
, ecf => [@ecfList]
, pimin => [@ebandlo]
, pimax => [@ebandhi]
)
or return exception();
doCommand(
'addattribute', set => $boxList
, attributename => 'BKGDSCRN'
, attributetype => 'boolean'
, booleanvalue => $flareScreened ? 'T' : 'F'
, attributecomment =>
"\"\\\"Was background screening of the eventlist applied?\\\"\""
)
or return exception();
&addAllExpIdKwds($boxList, \%expidKwds) or return exception();
# Step 5: Create spline background map for each image ----------------------
#
foreach my $i (0 .. $#rawImList)
{
my $bkgFile = $bkgImList[$i];
# info("DEBUG: maskImLongList [ $i ]".$maskImLongList[$i]);
doCommand(
'esplinemap'
, boxlistset => $boxList
, idband => $bandLongList[$i]
, imageset => $filtImList[$i]
, detmaskset => $maskImLongList[$i]
, withexpimage => 'true'
, expimageset => $flatExpImList[$i]
, withdetmask => 'true'
, withexpimage2 => 'false'
, fitmethod => 'spline'
# , nsplinenodes => 12
, nsplinenodes => 13
# , excesssigma => 2.5
, excesssigma => 3.0
, nfitrun => 4
# , withootset => 'true'
, withootset => ($instList[$i] eq 'epn') ? 'true' : 'false'
# , mlmin => 1.0
, scut => 0.002
, ooteventset => $imEventList[$i]
, bkgimageset => $bkgFile
, pimin => $ebandlo[$i]
, pimax => $ebandhi[$i]
#### last 2 aren't needed if $imEventList[$i] is filtered on the energy band.
)
or return exception();
doCommand(
'addattribute', set => $bkgFile
, attributename => 'BKGDSCRN'
, attributetype => 'boolean'
, booleanvalue => $flareScreened ? 'T' : 'F'
, attributecomment =>
"\"\\\"Was background screening of the eventlist applied?\\\"\""
)
or return exception();
&addExpIdKwds($bkgFile, \%expidKwds, $instList[$i]) or return exception();
} # end loop over found merged images
# Now make band 8 background maps. These are made by summing, for each instrument, the background maps for bands 1 to 5. NOTE that, formally speaking, there is a possibility that not all the 5 constituents may be summed, in which case the background map cannot truly be said to be for band 8. There's also nothing but convention to say that the bands don't overlap. So from a software engineering point of view this is not a very watertight bit of code. But we have to leave something to future generations!
#
# Band 9 (xid) bkg maps now also made here, and the lists of xid images, exp maps and bkg maps loaded.
my ( @xidRawImList, @xidExpImList, @xidBkgImList, $bkgFile );
foreach my $inst (qw(epn emos1 emos2))
{
next unless exists $band8constituents{$inst};
if (@{$band8constituents{$inst}} != 5)
{
info("There are not (as there should be) 5 constituents to the 'band 8' background map for instrument $inst.");
}
next unless(@{$band8constituents{$inst}});
$bkgFile = newFile(
class => 'intermediate'
, instrument => $inst
, band => 8
, content => 'merged background map'
, format => 'FITS'
);
my $tempImage = newFile(
class => 'intermediate'
, instrument => 'epic'
, content => 'temporaryImage'
);
doCommand(
'imweightadd', imagesets => $band8constituents{$inst}
, withweights => 'no'
, outimageset => $bkgFile
, tempset => $tempImage
)
or return exception();
my @bkgmaplistdtl = ( class => 'intermediate'
, content => 'background map list'
, instrument => $inst
, band => 8
, format => 'ASCII'
);
unless ( &ModuleUtil::keepSafe ( file_identifier => \@bkgmaplistdtl , text => $band8constituents{$inst} ) )
{
info('Unable to store band 8 image list for later use.');
return exception();
}
# Now for band 9 (xid):
my @xidBandsThisInst = @{$xidBands->{$inst}};
next unless exists $band9constituents{$inst};
if (@{$band9constituents{$inst}} != @xidBandsThisInst)
{
info("Incorrect number of constituents found for the XID background map for instrument $inst.");
}
next unless(@{$band9constituents{$inst}});
$bkgFile = newFile(
class => 'intermediate'
, instrument => $inst
, band => 9
, content => 'merged background map'
);
doCommand(
'imweightadd', imagesets => $band9constituents{$inst}
, withweights => 'no'
, outimageset => $bkgFile
, tempset => $tempImage
) or return exception();
push( @xidBkgImList, $bkgFile );
my $rawim = findFile(
class => 'intermediate'
, instrument => $inst
, band => 9
, content => 'merged image'
);
my $expim = findFile(
class => 'intermediate'
, instrument => $inst
, band => 9
, content => 'merged exposure map'
);
if ( $rawim && $expim )
{
push( @xidRawImList, $rawim );
push( @xidExpImList, $expim );
}
else
{
info("Cannot find all required images for instrument $inst, XID energy band.");
return exception();
}
} # end loop over instruments
# Step 6: Re-run eboxdetect with a background map -------------------------------------
#
my $mboxList = newFile(
class => 'product'
, instrument => thisInstrument
, content => 'EPIC observation box-map source list'
);
doCommand(
'eboxdetect'
, boxsize => 5
, imagesets => [@filtImList]
, expimagesets => [@expImList]
, bkgimagesets => [@bkgImList]
, detmasksets => [@band1MaskList]
, boxlistset => $mboxList
, likemin => 5
, usemap => 'true'
, withdetmask => 'true'
, withexpimage => 'true'
, nruns => 1
, ecf => [@ecfList]
, pimin => [@ebandlo]
, pimax => [@ebandhi]
)
or return exception();
doCommand(
'addattribute', set => $mboxList
, attributename => 'BKGDSCRN'
, attributetype => 'boolean'
, booleanvalue => $flareScreened ? 'T' : 'F'
, attributecomment =>
"\"\\\"Was background screening of the eventlist applied\\\"\""
)
or return exception();
&addAllExpIdKwds($mboxList, \%expidKwds) or return exception();
if ( fileExists( file => $mboxList )
&& numberFITSRows( file => $mboxList, extension => 'SRCLIST' )
> 0 )
{
my (@commandLine, %hrmatch);
# Step 7: Run emldetect on standard bands -----------------------------------
#
# EC: vignetting factors output by emldetect are currently
# wrong. The temporary fix is to run ep_updatevig to hack the
# correct columns in the table. Therefore we send the output
# of emldetect to intermediate/ instead of product/, and then
# place the hacked version output by ep_updatevig in product/.
# use this one temporarily for the emldetect call: remove when
# emldetect is updated in the future.
my $tempmlList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, content => 'EPIC observation ml source list'
);
# we still create this one for the final/corrected product
my $mlList = newFile(
class => 'product'
, instrument => thisInstrument
, content => 'EPIC observation ml source list'
);
my @srcMapList = @rawImList; ## added for psfdev
@srcMapList = map { s/merged/srcmap/i; $_; } @srcMapList; ## added for psfdev
@commandLine = (
imagesets => [@rawImList]
, expimagesets => [@expImList]
, bkgimagesets => [@bkgImList]
, sourceimagesets => [@srcMapList] ## added for psfdev
, boxlistset => $mboxList
, ecut => 15.0
, nmulsou => 2
, ecf => [@ecfList]
, determineerrors => 'true'
, mlmin => 6
, pimin => [@ebandlo]
, pimax => [@ebandhi]
# switch back to $mlList when emldetect is updated
# , mllistset => $mlList
, mllistset => $tempmlList
, dmlextmin => 4
, fitextent => 'true'
, withsourcemap => 'false' ## true for psfdev, false for sasdev
### GL has true, but this has no effect on source detection.
, psfmodel => 'ellbeta' ## ellbeta for psfdev, not present for sasdev
, withdetmask => 'true'
# , detmasksets => [@maskImList]
, detmasksets => [@band1MaskList]
, extentmodel => 'beta'
, withthreshold => 'true'
, threshcolumn => 'LIKE'
, threshold => 10.0
, withtwostage => 'true'
);
# Add hardness details to commandLine based on available instruments
%hrmatch = (
'emos1' => 'hrm1def'
, 'emos2' => 'hrm2def'
, 'epn' => 'hrpndef'
);
foreach my $inst (qw(epn emos1 emos2))
{
push( @commandLine
, ( $hrmatch{"$inst"}, [ @{ hardnessBands($inst) } ] )
);
}
doCommand( 'emldetect', @commandLine ) or return exception();
# If emldetect didn't find any sources, the table has no
# rows
my $nsrc = numberFITSRows(file => $tempmlList,
extension => "SRCLIST");
if ( $nsrc > 0 ) {
# Now fix up the vignetting columns, outputting the corrected version
# to $mlList. Remove this call once emldetect is updated.
doCommand( 'ep_updatevig',
useemllistset => 'yes',
pimin => '"200 500 1000 2000 4500 200 500 1000 2000 4500 200 500 1000 2000 4500"',
pimax => '"500 1000 2000 4500 12000 500 1000 2000 4500 12000 500 1000 2000 4500 12000"',
emllistset => $tempmlList,
emloutset => $mlList );
} else {
# If emldetect didn't find any sources, skip ep_updatevig and just copy to the final list
info("emldetect didn't seem to find any sources. Skipping ep_updatevig.");
copyFile(
source => $tempmlList
, destination => $mlList
);
}
&addAllExpIdKwds($mlList, \%expidKwds) or return exception();
# -----
my $backupFile1 = newFile(
class => 'intermediate'
, instrument => thisInstrument
, content => 'ML sourcelist dpssflag check 1'
);
copyFile(
source => $mlList
, destination => $backupFile1
);
# -----
# Now run emldetect on xid images:
if (@xidRawImList != @xidEcfList)
{
info("Wrong number of input files for XID run of emldetect");
}
else
{
my $mlListXid = newFile(
class => 'intermediate'
, instrument => thisInstrument
, content => 'ml XID source list'
);
@commandLine = (
imagesets => [@xidRawImList]
, expimagesets => [@xidExpImList]
, bkgimagesets => [@xidBkgImList]
, boxlistset => $mlList
, ecut => 15.0
, nmulsou => 1
, ecf => [@xidEcfList]
## GL says is obsolete parameter
, determineerrors => 'true'
, mlmin => 0.0
, pimin => [@xidEbandlo]
, pimax => [@xidEbandhi]
, mllistset => $mlListXid
, fitextent => 'no'
, fitposition => 'no'
, withsourcemap => 'false'
, psfmodel => 'ellbeta' ## ellbeta for psfdev, not present for sasdev
, withdetmask => 'true'
, detmasksets => [@band1MaskList]
, extentmodel => 'beta'
, xidfixed => 'yes'
, withthreshold => 'true'
, threshcolumn => 'LIKE'
, threshold => 10.0
);
# Add hardness and xid band details to commandLine based on available instruments
%hrmatch = (
'emos1' => 'hrm1def'
, 'emos2' => 'hrm2def'
, 'epn' => 'hrpndef'
);
my %xidmatch = (
'emos1' => 'xidm1def'
, 'emos2' => 'xidm2def'
, 'epn' => 'xidpndef'
);
foreach my $inst (qw(epn emos1 emos2))
{
push( @commandLine
, ( $hrmatch{"$inst"}, [ @{ hardnessBands($inst) } ] )
);
push( @commandLine
, ( $xidmatch{"$inst"}, [ @{ $xidBands->{$inst} } ] ) );
}
doCommand( 'emldetect', @commandLine ) or return exception();
doCommand(
'addattribute', set => $mlList
, attributename => 'BKGDSCRN'
, attributetype => 'boolean'
, booleanvalue => $flareScreened ? 'T' : 'F'
, attributecomment =>
"\"\\\"Was background screening of the eventlist applied?\\\"\""
)
or return exception();
# -----
my $backupFile2 = newFile(
class => 'intermediate'
, instrument => thisInstrument
, content => 'ML sourcelist dpssflag check 2'
);
copyFile(
source => $mlList
, destination => $backupFile2
);
# -----
}
}
# Create sensitivity map for band 8.
my $senseFile = newFile(
class => 'product'
, instrument => 'epic'
, band => 8
, content => 'EPIC OBSERVATION SENSITIVITY MAP'
);
doCommand(
'esensmap'
, expimagesets => [@expImList]
, bkgimagesets => [@bkgImList]
# , detmasksets => [@maskImList]
, detmasksets => [@band1MaskList]
, sensimageset => $senseFile
, mlmin => 6
###### Should be the same as for emldetect?? Or eboxdetect --likemin????
)
or return exception();
doCommand(
'addattribute', set => $senseFile
, attributename => 'BKGDSCRN'
, attributetype => 'boolean'
, booleanvalue => $flareScreened ? 'T' : 'F'
, attributecomment =>
"\"\\\"Was background screening of the eventlist applied?\\\"\""
)
or return exception();
# Generate OOTE combined mask ----------------------------------------
my @ootemaskarray = ();
my @minstlist = ();
my $instnum = 0;
foreach my $inst (qw(epn emos1 emos2)) {
$instnum++; # numerical code for instrument
my $instootemask = findFile(
class => 'intermediate'
, instrument => $inst
, band => 8
, format => 'FITS'
, content => 'Pileup mask'
);
if ( fileExists( file => $instootemask ) ) {
push( @ootemaskarray, $instootemask );
push( @minstlist, $instnum );
}
}
my $mlList = findFile(
class => 'product'
, instrument => thisInstrument
, content => 'EPIC observation ml source list'
);
# Check for existence of ML source list and run OOTE flagging
if ($mlList) {
# Sum number of OOTE masks
my $numootemasks = scalar @ootemaskarray;
my $pileupMaskSum = newFile(
class => 'intermediate'
, instrument => 'all'
, content => 'Pileup mask sum'
, band => 8
, format => 'FITS'
);
if ($numootemasks eq 0) {
info("No OOTE pileup masks found");
# return success();
} elsif ($numootemasks eq 1) {
info("One OOTE pileup mask found");
$pileupMaskSum = $ootemaskarray[0];
} elsif ($numootemasks eq 2) {
info("2 OOTE pileup masks found");
doCommand(
'ftimgcalc'
, outfile => $pileupMaskSum
, expr => '(A + B < 2. ? 0 : 1)'
, a => "A=$ootemaskarray[0]"
, b => "B=$ootemaskarray[1]"
);
} elsif ($numootemasks eq 3) {
info("3 OOTE pileup masks found");
doCommand(
'ftimgcalc'
, outfile => $pileupMaskSum
, expr => '(A + B + C < 3. ? 0 : 1)'
, a => "A=$ootemaskarray[0]"
, b => "B=$ootemaskarray[1]"
, c => "C=$ootemaskarray[2]"
);
}
# Run dpssflag on output merged source list
if ($numootemasks > 0) {
push( @ootemaskarray, $pileupMaskSum );
push( @minstlist, 0 );
doCommand(
'dpssflag'
, srclisttable => "$mlList:SRCLIST"
, extlimit => 1.5
, mfraclimit => 0.5
, brcolumn => 'RATE'
, brfactor => 65
, op1 => 'sqrt'
, op2 => 'mul'
, minrad => 10
, maxrad => 400
, extfactor => 3
, maxextrad => 200
, brlimit => 1000
, fixext => 160
, brexfraction => 0.4
, detmlfraction => 0.9
, badcolumn => 'CUTRAD'
, badfactor => 1
, withflag10 => 'yes'
, ootrgamasklist => [@ootemaskarray]
, minstlist => [@minstlist]
);
} else {
doCommand(
'dpssflag'
, srclisttable => "$mlList:SRCLIST"
, extlimit => 1.5
, mfraclimit => 0.5
, brcolumn => 'RATE'
, brfactor => 65
, op1 => 'sqrt'
, op2 => 'mul'
, minrad => 10
, maxrad => 400
, extfactor => 3
, maxextrad => 200
, brlimit => 1000
, fixext => 160
, brexfraction => 0.4
, detmlfraction => 0.9
, badcolumn => 'CUTRAD'
, badfactor => 1
, withflag10 => 'no'
);
}
}
# Raise success flag
return success();
}
sub addAllExpIdKwds
{
my $status = 1;
my ($fileName, $expidKwdsRef) = @_;
foreach my $inst (qw(epn emos1 emos2))
{
$status = &addExpIdKwds($fileName, $expidKwdsRef, $inst);
return $status if (!$status);
}
return $status;
}
sub addExpIdKwds
{
my $status = 1;
my ($fileName, $expidKwdsRef, $inst) = @_;
my %expidKwds = %{$expidKwdsRef};
my $prefix = 'PN'; # default
if ($inst eq 'emos1')
{
$prefix = 'M1';
}
elsif ($inst eq 'emos2')
{
$prefix = 'M2';
}
if (defined($expidKwds{$inst}))
{
my @kwd_list = @{$expidKwds{$inst}};
foreach my $i (0 .. $#kwd_list)
{
my $exp_id = $kwd_list[$i];
my $kwd_str = $prefix.'EXID'.sprintf("%2.2d", $i+1);
doCommand(
'addattribute', set => $fileName
, attributename => $kwd_str
, attributetype => 'string'
, stringvalue => $exp_id
)
or return exception();
}
}
return $status;
}
sub isInList
{
my ($element, @list) = @_;
my $status = 0;
return $status if (!defined($list[0]) || !defined($element));
foreach my $trial (@list)
{
if ($trial eq $element)
{
$status = 1;
last;
}
}
return $status;
}
1;