package EPICSourceProducts;
use strict;
use English;
use Carp;
use List::Util qw(min max);
use Regexp::Common qw(boolean);
use vars
qw(@ISA $VERSION $name $author $date $version @instList $numberOfStreams);
@ISA = qw(Module);
use ModuleResources;
use ModuleUtil;
use Astro::FITS::CFITSIO qw(:longnames :constants);
use FITSException;
use Background; # provides select_bkg() (and select_bkg_ebkgreg() for PN)
use Data::Dumper;
# Declare identity, version, author, date, etc.
$name = __PACKAGE__;
$VERSION = '1.74';
$version = $VERSION;
$author = 'Anja Schroeder,Ian Stewart,Duncan John Fyfe,Grant Denkinson,Duncan Law-Green,Ed Chapin, Jose V Perea';
$date = '2017-06-20';
#
# ChangeLog
# =========
#
# Version 1.74 - 2017-06-20 JVP
# ------------
#
# + Background region files for ANNULUS in MOS have now the same COMPONENT value (COMP=1) for all the region entries:
# For the ANNULUS and for the contaminating sources !CIRCLE
#
# Version 1.73 - 2017-03-09 JVP
# ------------
#
# + Contaminating source (red circles) Source ID labels were wrong in background region files (loop order instead the SourceID)
# + 'EPIC summary source list' TSERIES column is True only if any epic Timeseries exists
#
# Version 1.72 - 2017-01-30 JVP
# ------------
#
# + Small bug correction: when 'implotregion' task fails a log meesage is produced, but not Exception
#
# Version 1.71 - 2017-01-23 JVP
# ------------
#
# + Small bug correction in the index name of the last fork() loop: $j -> $i
# $BKG[$j] eq 'NONE' --> $BKG[$i] eq 'NONE'
#
# Version 1.70 - 2016-09-19 JVP
# ------------
#
# + Include PNG for 'EPIC SOURCE SPECTRUM'. PNG contains the thumbnail image of the source and background extraction regions.
# 'implotregions' task creates the thumbnail images.
#
# Version 1.69 - 2016-07-21 JVP
# ------------
#
# + SAS task for background region estimation 'ebkgreg' has been extended to MOS observations.
# + Module Background.pm modified
# + Source position is S/N optimizated (eregionanalyse) only for BKG_SHAPE = CIRCLE and skipped for ANNULUS.
#
# Version 1.68 - 2016-04-05 JVP
# ------------
#
# + OFFAX limited to maximum value from LOOKUP.FITS file (PPR#7333)
# Module 'List::Util qw(min max)' included to find min and max from array
# + Use fixed center in eregionanalyse. Source position from PSF fit instead 'eregionanalyse' (PCR#7307)
# + Skip epiclccorr when EXPOSURE keyword is missing also in background 'uncorrected source time series'
# + EPIC-pn : Do not ignore ~ 8 KeV spectral region (7.875:8.230) for spectrum products (PCR#7353)
# + SPECTRA & TSERIES flags: No Exception if there are not source spectra from any EPIC camera
#
# Version 1.67 - 2016-03-16 JVP
# ------------
#
# + Intermediate spectra with EXPOSURE absent or set to zero should not be copied to final products. PPR #7345
#
# Version 1.66 - 2016-03-03 JVP
# ------------
#
# + Flags columns SPECTRA & TSERIES = True in 'EPIC summary source list' for sources for which
# products are produced from at least 1 EPIC camera (epn,emo1/2). Final corrections.
#
# Version 1.65 - 2016-01-12 JVP
# ------------
#
# + Values captured from the stdout in 'ecoordconv' and 'eregionanalyse' can be now non-integer
# + Flags columns SPECTRA & TSERIES = True in 'EPIC summary source list' only for @selectedSources
# + EPIC PN: If background is not possible for any source (ex. rb < 3 pix)
# => no products => skip to next source
# + Flags columns SPECTRA & TSERIES = True in 'EPIC summary source list' for sources for which
# products are produced from at least 1 EPIC camera (epn,emo1/2)
#
# Version 1.64 - 2015-04-10 JVP
# ------------
#
# + New SAS task for background region searching in PN observations: 'ebkgreg'.
# New subroutine called : 'select_bkg_ebkgreg'. It is defined in Background.pm module.
# Background region obtained in the same CCD where the source is located.
# MOS observations keep the previous method (iterative tests arround the source).
#
# Version 1.63 - 2015-04-10 JVP
# ------------
#
# + Skip epiclccorr when EXPOSURE keyword is missing in 'uncorrected source time series'. PPR #7281
#
# Version 1.62 - 2014-11-03 JVP
# ------------
#
# + Units of source extraction radius in output spectral products corrected (arsec -> 0.05 arcsec)
#
# Version 1.61 - 2014-07-03 JVP
# ------------
#
# + Iteration over x2, x4 binning factors in light curves generation removed (commented).
# They only were 'intermediate' products.
#
# Version 1.60 - 2014-06-25 JVP
# ------------
#
# + Call 'grppha' replaced by 'specgroup'. Bad points selected by energy range, instead channel ranges.
# FFT Timeseries plot (efftplot) produced from final product (filename label corrected) and only generated for $specTotalCounts > 100.
# Re-naming ANCRFILE, RESPFILE and BACKFILE keywords with only one 'addattribute' call.
#
# Version 1.59 - 2014-05-28 JVP
# ------------
#
# + Call to 'lccorr_pcms' replaced by 'epiclccorr' in lccorr_succeededN command.
#
# Version 1.58 - 2014-05-16 JVP
# ------------
#
# + Call to 'lccorr_pcms' replaced by 'epiclccorr' in lccorr_succeeded command.
# (pending to replace in lccorr_succeededN command).
#
# Version 1.57 - 2013-02-24 EC
# ------------
#
# + If especget fails to produce spectrum, immediately skip to next source
#
# Version 1.56 - 2013-12-09 EC
# ------------
#
# + Parallelize all of the main loops using forks
#
# Version 1.55 - 2013-06-11 EC
# ------------
#
# + Fix filename for light curves (add band=>8)
#
# Version 1.54 - 2013-06-11 EC
# ------------
#
# + Merge changes from SOC version.
#
# version 1.53 2012-11-06 DLG
# ------------
#
# - Correct ANCRFILE keyword to point to ARF product file (.FTZ)
#
# version 1.52 2012-10-25 DLG
# ------------
#
# - Generate source ARF as intermediate, copy to product for >100 counts
# - Move generation of ASCII region files inside counts>100 loop
#
# version 1.51 2012-10-01 DLG
# ------------
#
# - Revert 1.50 for SSP spectra.
#
# version 1.50 2012-09-26 DLG
# ------------
#
# + Add FITS source and background region extensions to SSP spectra and timeseries
#
# version 1.49 2012-09-11 DLG
# ------------
#
# + Change threshold illumination factor to 0.70
#
# version 1.48 2012-09-11 DLG
# ------------
#
# + Fix origin filename for timeseries PDF plot (product, not intermediate)
#
# version 1.47 2012-08-20 DLG
# ------------
#
# + set SSP offset for SW mode from 70 to 45 (larger radial steps to speed processing)
#
# version 1.46 2012-08-14 DLG
# ------------
#
# + replace flare GTI with optimised flare GTI if found
#
# version 1.45 2012-08-02 DLG
# ------------
#
# + write source extraction radius to output spectral products
#
# version 1.44 2012-07-26 DLG
# ------------
#
# + clamp eregionanalyse minimum radius to 240 (12")
#
# version 1.43 2012-07-04 DLG
# ------------
#
# + eregionanalyze check for zero source counts in extraction region
#
# version 1.42 2012-07-04 DLG
# ------------
#
# + Test lowering global illumination threshold to 0.85
#
# version 1.41 2012-07-03 DLG
# ------------
#
# + Generate intermediate lightcurves at x2 and x4 binning
#
# version 1.40 2012-06-27 DLG
# ------------
#
# + call ekstest twice to generate chisq and fvar values
#
# version 1.39 2012-06-27 DLG
# ------------
#
# + Combine extraction detmask with RGA/OOT mask, where found
#
# version 1.38 2012-06-12 DLG
# ------------
#
# + Only create product bkg spectra, timeseries + associated plots where
# SSP spectrum has > 100 valid counts
#
# version 1.37 2012-04-30 DLG
# ------------
#
# + Propagate source flags to source specific products (new code)
#
# version 1.36 2012-04-30 DLG
# ------------
#
# + Add ekstest fracvartest => 'yes'
#
# version 1.35 2012-03-30 DLG
# ------------
#
# + skip spectral product generation on <100 good counts
#
# version 1.34 2012-03-27 DLG
# ------------
#
# + Generate intermediate file with total good spectral counts (QUALITY=0)
#
# version 1.33 2012-03-26 DLG
# ------------
#
# + Check for faint source counts in main optimisation loop to avoid xspec errors
#
# version 1.32 2012-03-05 DLG
# ------------
#
# + Check for existence of ${cam}_CTS column in $sspSrcInput before filtering
#
# version 1.31 2012-02-13 DLG
# ------------
#
# + Modify ASCII DS9 region file output to reflect new optimised regions
# + Apply optimised extraction regions to lightcurve generation
#
# version 1.30 2011-10-31 DLG
# ------------
#
# + implementing background-optimised spectral extraction
#
# version 1.292 2013-02-03 EC
# ------------
#
# + Remember to add source attributes even if skipping grppha
#
# version 1.291 2013-01-31 EC
# ------------
#
# + Skip grppha for spectra missing EXPOSURE keyword after especget
#
# version 1.29 2009-10-16 GD
# ------------
#
# + Comment out changes to propagate flags to timeseries and
# spectral source specific products while awaiting esources SCR 286
#
# version 1.28 2009-09-09 GD
# ------------
#
# + Read automatic source flags ready to propagate to timeseries and
# spectral source specific products
#
# version 1.27 2009-09-03 DLG
# ------------
#
# + Manually set AVRATE, CHISQUAR, CHI2PROB to -1 and N_POINTS to 0
# in timeseries output prior to ekstest run.
#
#
# version 1.26 2009-01-26 DLG
# ------------
#
# + Corrected ANCRFILE and BACKFILE headers in output filename to reflect compressed FITS files
#
#
# version 1.25 2006-10-13 DJF
# ------------
#
# + Make sure each image exists before we try to access fits keywords,
# + Make sure we have some images to act on aswell.
#
# version 1.24 2006-09-20 DJF
# ------------
#
# + Record use of exposures for source specific product creation in the database.
#
# version 1.23 2006-09-15 DJF
# ------------
#
# + evaluateRules added test on ExpChooser. Exposures ignored by ExpChooser were being processed. This should not happen.
#
# version 1.22 - 2006-09-14 DJF
# ------------
#
# + Construct short unique filename to replace hardcoded xpsec.ps for xpsec postscript output.
# + Apply isEmptyImage test to band 1-5 images. If any are empty do not proceed with source specific product creation.
#
# version 1.21 - 2006-09-05 DJF
# ------------
#
# + Use fmodhead to re-write ANCRFILE keyword in spectrum file SPECTRUM extension. the extension .FIT becomes .FTS to maintain post-compression usability.
#
# version 1.20 - 2006-08-08 ACS
# ------------
#
# + changed line width in xspec plot from 3 to 6 (for better conversion into gif-files)
#
# version 1.19 - 2006-07-13 DJF
# ------------
#
# + Removed # from source spectrum title - xpsec doesn't like it.
#
# version 1.18 - 2006-07-13 DJF
# ------------
#
# + Added decimal source number and data source filename excluding leading directory to title of spectral plots
#
# version 1.17 - 2006-06-01 ACS
# ------------
#
# + Removed filename from label in spectrum plot
#
# version 1.16 - 2006-04-11 DJF
# ------------
#
# + Changed use of echeckregion. We cannot use the return code
# so the task must write output to a file and we must parse that file.
#
# version 1.15 - 2006-03-31 DJF
# ------------
#
# + Modified search for canned rmf. Fail if one is not found.
# + Use of TLMIN/TLMAX can leave light curves with all zero 'RATE'. Now check for this and avoid doing anything further with them.
# + EXPOSU?? sort was looking for CCD 1 .. 12 for both MOS and PN. Now restricted to 1 .. 7 for MOS.
# + Test of BKGDSCRN for flare background screening was testing if keyword was set or not, not if it were T or F (or Y/N ...)
#
# version 1.14 - 2006-03-15 DJF
# ------------
#
# + Error check in ML_ID_SRC selection called too early. Tis stopped the test completing. Moved.
#
# version 1.13 - 2006-03-15 DJF
# ------------
#
# + Removed attempt to make SSP for PN small window mode. This needs more careful consideration. Without PN SW mode source detection it can be impossible to match SSP selected sources with emldetect source list entries.
#
# version 1.12 - 2006-03-10 DJF
# ------------
#
# + fmodhead template creation moved out of if(){} statement scope
# so it can be used for additional keyword setting.
#
# version 1.11 - 2006-03-09 ACS
# ------------
#
# + Changed writing of source-specific keywords into the header of the spectrum and arf-files
# from 'addattribute' to 'fmodhead'; (also corrected keyword comments for addattribute in @kwdCommentList)
#
# version 1.10 - 2006-03-03 ACS
# ------------
#
# + Added very crude solution to extract source-level products for PN SW when the source lists
# have no PN entries at all (and esources therefore indicates that no PN sources
# will be extracted). Now _all_ sources in the SSP source list will be used,
# and echeckregion will through out all sources _not_ in the detector window.
#
# version 1.09 - 2006-03-02 ACS
# ------------
#
# + Removed sections about TSART and TSTOP keywords (see v 1.07)
# since it was done in the wrong place with the wrong values;
# moved this part to EPICSources.pm
#
# version 1.08 - 2006-03-02 ACS
# ------------
#
# + Add echeckregion after making source extraction region,
# will skip the source if it lies outside the detector window
#
# version 1.07 - 2006-02-21 DJF
# ------------
#
# + Write TSART and TSTOP keywords from esources output to
# lccorr output as TLMIN1 and TLMAX1 ahead of calling evselect.
# This is to cause timeseries to be aligned on these values.
#
#
# version 1.06 - 2006-02-20 DJF
# ------------
#
# + Sources were being selected with an foo() and info() if X
# This depends on info() returning true which is not guarenteed.
# changed to a more normal if(X) { foo(); info(); }
#
# version 1.05 - 2006-02-07 DJF
# ------------
#
# + Sort events for MOS and PN. There are a few ODF where MOS frames are not in order.
# This causes lccorr problems so we sort them first.
#
# version 1.04 - 2006-01-25 DJF
# ------------
#
# + Fixed unnecessary (and problematic) quoting of ftool parameters.
#
# version 1.03 - 2005-12-06 DJF
# ------------
# + Raise an exception if catASCIIFiles fails.
#
#
# version 1.02 - 2005-11-30 IMS
# ------------
# + Added --maketimecolumn=yes and --makeratecolumn=yes to the new evselect calls.
#
# version 1.01 - 2005-11-30 IMS
# ------------
# + Changed single quotes to double in --expression values for the two new evselect calls.
#
# version 1.0 - 2005-11-29 IMS
# ------------
# + I've replaced etimeget by 2 calls of evselect. This is because I find that etimeget does filtering on FLAG and PATTERN without being told to. We need to keep control of this within the pipeline.
# + SSPs now not made unless the exposure was flare-screened (requested by ACS 20051129).
# + hasFITSKeyword added to the list of used subroutines of ModuleResources.
#
# version 0.24 - 2005-11-25 IMS
# ------------
# + --updateexposure=no added to 2nd evselect call.
#
# version 0.23 - 2005-11-25 IMS
# ------------
# + Existence of spectrum PS file now checked before PStoPDF is invoked.
# + I'd forgotten to copy the sorted event list back to $forTimeFilteredEvList. D'oh!!
# + The DS9 region product is now not made until after the time series and the spectrum. There's no point in making a DS9 region if there is neither useful spectral or TS data.
#
# version 0.22 - 2005-11-24 IMS
# ------------
# + TS products now labelled as band 8.
#
# version 0.21 - 2005-11-23 IMS
# ------------
#
# + Replaced the first evselect call (filters the ml source list) by fselect. This is because evselect failed with error 'Accessing data SRCLIST.CCDNR of type Int32 as int8 is not supported'.
#
# version 0.20 - 2005-11-23 IMS
# ------------
#
# + Altered it so that the emldetect source list (suitably filtered) is supplied to region instead of the SSP-filtered list (made within EPICSources from the srcmatch/eposcorr/evalcorr output list). This requires some careful accounting to make sure that we are always talking about the same source.
# + Now tests for existence of an EXPOSUnn extension before attempting to sort it (see ChangeLog entry 0.18).
# + Added hasFITSExtension to the list of subroutines called from ModuleResources.
#
# version 0.19 - 2005-11-22 IMS
# ------------
#
# + Changed the order of the ds9 components in catASCIIFiles so the white circles don't get overdrawn with cyan.
#
# version 0.18 - 2005-11-21 IMS
# ------------
#
# + Module now tests that time series have more than 0 non-null RATE values (efftplot etc not called if not).
# + EXPOSU extensions of pn event lists are now sorted.
#
# version 0.17 - 2005-11-16 IMS
# ------------
#
# + Added a couple more infos.
#
# version 0.16 - 2005-11-16 IMS
# ------------
#
# + Module no longer fails if efftplot fails - just the subsequent psToPdf is not called.
#
# version 0.15 - 2005-11-15 IMS
# ------------
#
# + Added a print of the rmf file name as read from the temp spectrum, in an attempt to track down why the module can't work out the path name.
# + Added carriage returns to the fmodhead template file!
#
# version 0.14 - 2005-11-14 IMS
# ------------
#
# + Incorporating slight changed DJF made to intermediate file names.
# + Somewhat more care is now taken in juggling the ARF, RMF and bkgSpectrum file vs path names in the source-spectrum SPECTRUM extension header. The point is that xspec needs the path name; but after the xspec call, because pcms path names are meaningless to the user, just the bare file names are substituted, using now a second grppha call.
#
# version 0.13 - 2005-11-03 IMS
# ------------
#
# + Label of spectrum pdf plot set to 'source minus background' instead of the default.
# + fmodhead template file (introduced in v0.9) was left with just place-holder values! These have now been corrected.
#
# version 0.12 - 2005-10-31 IMS
# ------------
#
# + Reversed extraction of bare file name from the path name for the rmf file. I think especget writes this bare name to the header already; don't have to fiddle with it. But check this.
# + Changed content string for source ds9 region product to be consistent with FileMap.
#
# version 0.11 - 2005-10-28 IMS
# ------------
#
# + Modified the xspec command file a little in order to change the y-axis scale of spectral plots from linear to log (request from ACS).
# + Changed the grppha invocation such that the bare file names rather than path names of the arf, rmf and bkg spectrum files are written as keywords to the src spectrum headers.
#
# version 0.10 - 2005-10-27 IMS
# ------------
#
# + Now setting the values of lccorr_pcms parameters srcra and srcdec; choice for --srcweightstyle now changed from 'events' to 'deduce'; new parameter --srcdeducestyle='events' now set.
#
# version 0.09 - 2005-10-14 IMS
# ------------
#
# + Shortened the content string for the SSP source list.
# + Call to lccorr_pcms now uses --srcweightstyle='events' rather than the default 'psf'. It is hoped that this will cure lccorr_pcms's remaining seg fault problems.
# + Added numberFITSRows, readFITSKeyword to list of used subroutines of ModuleResources.
# + elcplot and efftplot are now run on $correctedRateSet (the final product TS) instead of $bkgSubtractedRateSet, which is an intermediate file.
# + The first call to addattribute (the one which writes to the time series) has temporarily been replaced by a call to fmodhead. This is because, whereas elcplot (correctly) expects these keywords to be in the RATE extension header, addattribute is not yet capable of writing them to anywhere but the primary header. (Note that this ought not to stuff up efftplot since this task looks in both places for the kwds.)
# + Added new product (ds9 source regions).
# + Replaced a few instances of $$srcIdCol[$i] by $id.
#
# version 0.08 - 2005-09-23 RGW
# ------------
#
# + reinstate chi-square test in ekstest
#
#
# version 0.07 - 2005-09-23 IMS
# ------------
#
# + Various kwds get written to the SSP files. However they were all being written to :PRIMARY. These calls to addattribute have now been changed so that the kwds are written to the header of the main binary-table extension (= RATE for light curves, SPECTRUM for spectra and SPECRESP for the arf).
# + The RA_CORR and DEC_CORR columns are now being read from the $sspSrcFile. Previously the RA and DEC values were substituted.
# + Call to efftplot was commented out (?) - now restored.
#
# version 0.06 - 2005-09-14 IMS
# ------------
#
# + Added parameters --tempset and --tempeventset to lccorr_pcms call.
# + Failure in lccorr_pcms now does not cause the module to fail.
#
# version 0.05 - 2005-09-14 IMS
# ------------
#
# + Added extra parameters --srcregiontab and --bkgregiontab to the call to lccorr_pcms.
#
# version 0.04 - 2005-09-09 RGW
# ------------
#
# + fixed arguments to addattribute (integer, not int)
# + specify GTI block name to gtialign
# + use pre-defined content name for time series & plots
# + work around short filename limits in XSPEC
# + grppha requires "!" on output file to force overwrite
# + fixed parameter names on especget
# + not necessary to pass src number as a hex string to newFile
# + added a trap to test for lightcurves with zero bins (which can
# happen in observations which contain both full-frame and windowed
# exposures)
# + complete gracefully when no sources have been selected to product creation
# + added trap to ignore empty flare GTIs
#
# Version 0.03 - 2005-09-05 IMS
# ------------
#
# + Now calling lccorr_pcms, with 8 additional parameters recording the event selections, instead of lccorr.
#
# version 0.02 - 2005-08-24 RGW
# ------------
#
# + temporary hacks to get it running within the PCMS
#
# version 0.01 - 2005-06-09 IMS
# ------------
#
# + First draft. Anja's script converted to module format.
#
# Declare list of instruments this module is interested in
@instList = qw(epn emos1 emos2);
# Number of streams
sub numberOfStreams
{
return numberOfExposures();
}
# Rules method
sub evaluateRules
{
# If upstream module has been ignored, so if this one
return ignore() if ignored(module => 'ExpChooser'
,instrument => thisInstrument
,stream => thisStream
);
return ignore()
if ignored( module => 'EPICSources'
, instrument => 'all'
, stream => 1
);
return start()
if complete( module => 'EPICSources'
, instrument => 'all'
, stream => 1
);
}
# Action method
sub performAction
{
info("Module version number: $version");
my $exp_id = exposureID( instrument => thisInstrument
, stream => thisStream
);
info("Exposure ID: $exp_id");
my $flareScreened = undef;
# Make sure there are no empty band images and flare background screening has been applied.
my @imgcount;
foreach my $band ( 1 .. 5 )
{
my $image = findFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, band => $band
, content => 'EPIC image'
, 'format' => 'FITS'
);
next unless ($image && fileExists(file =>$image));
if (!defined($flareScreened))
{
if ( hasFITSKeyword( file => $image
, extension => 'PRIMARY'
, keyword => 'BKGDSCRN'
)) {
$flareScreened = readFITSKeyword(
file => $image
, extension => 'PRIMARY'
, keyword => 'BKGDSCRN'
);
$flareScreened = ($flareScreened =~ /$RE{boolean}{true}/i)
? 1
: 0
;
if (defined $flareScreened && !$flareScreened)
{
info("Exposure was not flare screened");
info("Source-specific products will not be made.");
return success();
}
}
}
my @fopt = ( class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, format => "ASCII"
, band => $band
, content => "Filtered image statistics"
);
my $stats = findFile( @fopt );
unless ($stats )
{
$stats = newFile( @fopt );
}
if ( ModuleUtil::isImageEmpty($image, $stats) )
{
info("$image is empty");
info("Excluding this exposure from source specific product creation.");
return success();
}
push @imgcount , $image;
}
unless (@imgcount)
{
info("This exposure has no images.");
info("Excluding this exposure from source specific product creation.");
return success();
}
# Get the ml source list:
#
my $mlList = findFile(
class => 'product'
, instrument => 'epic'
, content => 'EPIC observation ml source list'
);
return success() unless $mlList;
# get automatic source flags to be propagated - GD
# my $autoSrcFlags = readFITSColumn( file => $mlList
# , extension => 'SRCLIST'
# , column => 'FLAG'
# );
# my $autoSrcFlagsComment = readFITSComment( file => $mlList
# , extension => 'SRCLIST'
# , keyword => 'FLAG'
# );
# defined ($$autoSrcFlags[0]) && info("Automatic source flags: $$autoSrcFlags[0].");
# defined ($autoSrcFlagsComment) && info("Automatic source flags comment: $autoSrcFlagsComment.");
# my $autoSrcFlagsComment = "Automatic source flags";
# Filter the ml source list:
my $selExpression;
if (thisInstrument eq 'emos1')
{
$selExpression = "(DET_ML>15 && (ID_BAND==0)&&(ID_INST==2))";
}
elsif (thisInstrument eq 'emos2')
{
$selExpression = "(DET_ML>15 && (ID_BAND==0)&&(ID_INST==3))";
}
else # assume 'epn'
{
$selExpression = "(DET_ML>15 && (ID_BAND==0)&&(ID_INST==1))";
}
my $filteredMlList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, content => 'Filtered ML src list'
);
doCommand(
'fselect'
, infile => $mlList."[SRCLIST]"
, outfile => $filteredMlList
, expr => $selExpression
)
or return exception();
my $mlSrcIdColInOrig = readFITSColumn(
file => $filteredMlList
, extension => "SRCLIST"
, column => "ML_ID_SRC"
);
# Find the source-specific-products source list:
#
my $sspSrcInput = findFile(
class => 'intermediate'
, instrument => 'epic'
, content => 'SSP source list'
);
return success()
if ( !defined($sspSrcInput) || $sspSrcInput eq ''
|| !fileExists( file => $sspSrcInput )
);
# check for zero-length SSP source list
my $nsrc=numberFITSRows(file => $sspSrcInput,
extension => "SRCLIST");
info("SSP Source list contains $nsrc sources");
return success() if $nsrc==0;
# Strip out sources with null counts in this instrument
my $cam;
if (thisInstrument eq 'emos1') {
$cam = 'M1';
} elsif (thisInstrument eq 'emos2') {
$cam = 'M2';
} else {
$cam = 'PN';
};
# Generate name of source list counts column
my $cts_col = $cam."_CTS";
# Generate name of source list offaxis angle column
my $offax_col = $cam."_OFFAX";
# Filter SSP list for nonzero counts in camera
info("Filtering SSP list $sspSrcInput for $cts_col > 0");
my $sspSrcFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Filtered SSP source list'
, format => 'FITS'
) or return exception();
doCommand(
'fselect'
, infile => $sspSrcInput."[SRCLIST]"
, outfile => $sspSrcFile
, expr => "$cts_col > 0"
);
my $nsrc=numberFITSRows(file => $sspSrcFile,
extension => "SRCLIST");
info("Filtered Source list contains $nsrc sources");
return success() if $nsrc==0;
my $standardFlareGTI = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'flare GTI'
);
my $optFlareGTI = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Optimised flare GTI'
);
my $alignedGtis;
my $flareGtifile;
if (fileExists(file => $optFlareGTI)) {
$flareGtifile = $optFlareGTI;
} else {
$flareGtifile = $standardFlareGTI;
}
my $flareGtiExists = 0;
if (fileExists(file => $flareGtifile) &&
numberFITSRows(file => $flareGtifile,
extension => "STDGTI")) {
$flareGtiExists = 1;
}
my $forTimeFilteredEvList = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Event list for light curves'
);
if ( !defined($forTimeFilteredEvList)
|| $forTimeFilteredEvList eq ''
|| !fileExists( file => $forTimeFilteredEvList )
) {
info('No light curve event list. Nothing to do.');
return success();
}
### why is this 'success' but for the spectrum evlist 'exception'??
my $forSpecFilteredEvList = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Event list for spectra'
);
unless ( $forSpecFilteredEvList
&& fileExists( file => $forSpecFilteredEvList )
) {
info('No spectra event list. Nothing to do.');
return success();
}
my $flagBitKwdName = 'PN_P_BIT';
if (thisInstrument eq 'emos1')
{
$flagBitKwdName = 'M1_P_BIT';
} elsif (thisInstrument eq 'emos2')
{
$flagBitKwdName = 'M2_P_BIT';
}
my $flagBit = -1; # default
$flagBit = readFITSKeyword(
file => $sspSrcFile
, extension => "SRCLIST"
, keyword => $flagBitKwdName
);
return exception() if ($flagBit < 0);
my $flagMask = 1 << $flagBit;
# Set up the event selection parameter values for lccorr_pcms:
my ($patternHi, $eventFlagMask, $rawylo, @numCcds);
if (thisInstrument eq 'epn')
{
$patternHi = 4;
$eventFlagMask = '0xfffffef';
$rawylo = 0;
@numCcds = ( 1 .. 12 );
}
else # assume mos
{
$patternHi = 12;
$eventFlagMask = '0x766ba000';
$rawylo = 0;
@numCcds = ( 1 .. 7 );
}
# NOTE! The flag masks should be the same as the ones actually used in EPICSources.
my @gtiTabList;
foreach my $i (@numCcds)
{
push @gtiTabList, sprintf('%s:STDGTI%02d' , $forTimeFilteredEvList , $i);
}
# --------------------------------------------------------------------------------
my ($timeBinSizeCol, $srcIdCol, $srcIdHexCol, $mlSrcIdColInSSP);
my ($srcRaCol, $srcDecCol, $srcCorrRaCol, $srcCorrDecCol, $srcCtsCol, @srcOffaxCol);
my ($flagCol, $epFlagCol, $pnFlagCol, $m1FlagCol, $m2FlagCol);
my @selectedSources;
if ($sspSrcFile)
{
$timeBinSizeCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "TBIN_WIDTH"
);
$srcIdCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "SRC_NUM"
);
$srcIdHexCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "SRC_NUM_HEX"
);
$mlSrcIdColInSSP = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "ML_ID_SRC"
);
$srcRaCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "RA"
);
$srcDecCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "DEC"
);
$srcCorrRaCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "RA_CORR"
);
$srcCorrDecCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "DEC_CORR"
);
$flagCol = readFITSColumn(
file => $sspSrcFile
, extension => 'SRCLIST'
, column => 'PROCESS'
);
$epFlagCol = readFITSColumn(
file => $sspSrcFile
, extension => 'SRCLIST'
, column => 'EP_FLAG'
);
$pnFlagCol = readFITSColumn(
file => $sspSrcFile
, extension => 'SRCLIST'
,column => 'PN_FLAG'
);
$m1FlagCol = readFITSColumn(
file => $sspSrcFile
, extension => 'SRCLIST'
, column => 'M1_FLAG'
);
$m2FlagCol = readFITSColumn(
file => $sspSrcFile
, extension => 'SRCLIST'
, column => 'M2_FLAG'
);
$srcCtsCol = readFITSColumn(
file => $sspSrcFile
, extension => 'SRCLIST'
, column => "$cts_col"
);
# info("DEBUG: Got srcCtsCol[0]=".$$srcCtsCol[0]);
for ( my $i=0 ; defined $$flagCol[$i] ; $i++ )
{
return exception() if ( ! ( defined($$timeBinSizeCol[$i])
&& defined($$srcIdCol[$i])
&& defined($$srcIdHexCol[$i])
&& defined($$mlSrcIdColInSSP[$i])
));
if (($$flagCol[$i] & $flagMask) > 0)
{
push( @selectedSources, $i );
info("select source $i");
}
}
}
my %mlId_to_rowNum = ();
foreach my $i (@selectedSources)
{
#info("row check i=$i , \$mlSrcIdColInSSP->[$i] = $mlSrcIdColInSSP->[$i] , j = ( 0 .. $#$mlSrcIdColInOrig ) ");
foreach my $j (0 .. $#$mlSrcIdColInOrig)
{
#info("row check j=$j , \$mlSrcIdColInOrig->[$j] = $mlSrcIdColInOrig->[$j] ");
if ($mlSrcIdColInOrig->[$j] == $mlSrcIdColInSSP->[$i])
{
$mlId_to_rowNum{ $mlSrcIdColInSSP->[$i] } = $j+1;
}
}
unless ( exists( $mlId_to_rowNum{ $mlSrcIdColInSSP->[$i] } ) && defined( $mlId_to_rowNum{ $mlSrcIdColInSSP->[$i] } ) )
{
info("Can't find ML_ID_SRC value $mlSrcIdColInSSP->[$i] in $filteredMlList");
return exception();
}
}
# --------------------------------------------------------------------------------
# + DLG: add X,Y coords to source list
my ($mlIdCol, $mlRaCol, $mlDecCol, $mlCtsCol, $mlOffaxCol);
# get EPIC image
my $epicImage = findFile(
class => "product"
, instrument => thisInstrument
, exp_id => $exp_id
, band => 8
, content => 'EPIC image'
, format => 'FITS'
) or return exception();
my $summaryList = findFile(
class => "product"
, instrument => 'epic'
, content => 'EPIC summary source list'
, format => 'FITS'
) or return exception();
#-------------------------------------------------------------------------------------------------------------------
# Get SPECTRA and TSERIES flags columns from source list
my %src_spectra_tseries = readFITStoHash($summaryList);
# Flag Slected Sources to SPECTRA = TSERIES = TRUE in the %src_spectra_tseries hash
my @selectedSources_src_num; # @selectedSources_src_num are the SRC_NUM of the selectedSources
foreach my $i (@selectedSources){ # @selectedSources are the indexes of the SRC_NUM selectedSources
my $src = $$srcIdCol[$i];
push (@selectedSources_src_num, $src);
@selectedSources_src_num = sort { $a <=> $b } @selectedSources_src_num;
}
info("DEBUG for SPECTRA. selectedSources_src_num = @selectedSources_src_num ");
foreach my $source ( sort { $a <=> $b } keys(%src_spectra_tseries) ){
if ( $source ~~ @selectedSources_src_num ){ # If $id source is in @selectedSources list => True, True
@{$src_spectra_tseries{$source}} = (1,1); # ( SPECTRA, TSERIES ) = ( True, True )
#info("Source SRC_NUM = $source flagged to SPECTRA, TSERIES = True, True , because src is a selectedSources");
} else { # other sources NOT in @selectedSources => False, False
@{$src_spectra_tseries{$source}} = (0,0); # ( SPECTRA, TSERIES ) = ( False, False )
#info("Source SRC_NUM = $source flagged to SPECTRA, TSERIES = False, False , because src is NOT a selectedSources");
}
}
# # ==== DEBUG ==========================================================
# foreach my $source ( sort { $a <=> $b } keys(%src_spectra_tseries) ){
# info("DEBUG Reading SPECTRA hash after 1st re-writing: $source - @{$src_spectra_tseries{$source}}");
# }
# # ==== DEBUG ==========================================================
#-------------------------------------------------------------------------------------------------------------------
# Generate filtered source list
# remove sources with 0 or -ve counts
my $filteredSrcList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Column filtered source list'
, format => 'FITS'
);
doCommand(
'fselect'
, infile => $summaryList
, outfile => $filteredSrcList
, expr => "$cts_col > 0"
) or return exception();
# Get number of rows in filtered source list
my $num_sources = numberFITSRows(
file => $filteredSrcList
, extension => 'SRCLIST'
);
info("$filteredSrcList: found $num_sources with ${cam}_CTS > 0");
# Get number of sources in SSP list
my $num_ssp = numberFITSRows(
file => $sspSrcFile
, extension => 'SRCLIST'
);
# Generate column filtered source list
my $subSrcList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Column subset source list'
, format => 'FITS'
) or return exception();
my $columns = "SRC_NUM,$cts_col,RA,DEC,$offax_col";
doCommand(
'ftabcopy'
, infile => $filteredSrcList
, outfile => $subSrcList
, columns => $columns
# , rows => "1-$num_sources"
, rows => "-"
) or return exception();
# Read columns from filtered source list into arrays
$mlIdCol = readFITSColumn(
file => $subSrcList
, extension => "SRCLIST"
, column => 'SRC_NUM'
);
$mlRaCol = readFITSColumn(
file => $subSrcList
, extension => "SRCLIST"
, column => "RA"
);
$mlDecCol = readFITSColumn(
file => $subSrcList
, extension => "SRCLIST"
, column => "DEC"
);
$mlCtsCol = readFITSColumn(
file => $subSrcList
, extension => "SRCLIST"
, column => "$cts_col"
);
$mlOffaxCol = readFITSColumn(
file => $subSrcList
, extension => "SRCLIST"
, column => "$offax_col"
);
# Walk RA, DEC values, run coordinate conversion
my ($i, $j);
my (@allX_raw, @allY_raw, @all_CCD);
my $sccd;
for ($i = 0; $i < $num_sources; $i++) {
my $sc = $i+1;
my ($xcoord, $ycoord, $sccd) = coordConv($epicImage, $$mlRaCol[$i], $$mlDecCol[$i]);
info("coordConv returned - source $sc of $num_sources: x: $xcoord y: $ycoord CCD: $sccd");
# Build detector X,Y columns
push @allX_raw, int($xcoord);
push @allY_raw, int($ycoord);
push @all_CCD, int($sccd);
};
# Copy offaxis angles from ML source list to SSP source list
my %offax_hash;
@offax_hash{@$mlIdCol} = @$mlOffaxCol;
my @srcCCDCol;
for ( $j = 0; $j < $num_ssp; $j++ ) {
$srcOffaxCol[$j] = $offax_hash{$$srcIdCol[$j]};
$srcCCDCol[$j] = $all_CCD[$$srcIdCol[$j]];
};
# Write detector columns to source list
my $data;
$$data{X} = [@allX_raw];
$$data{Y} = [@allY_raw];
$$data{-column}{X} = { 'width' => 4,
'colnum' => 6,
'ttype' => 'X',
'repeat' => 1,
'typecode' => 42
};
$$data{-column}{Y} = { 'width' => 4,
'colnum' => 7,
'ttype' => 'Y',
'repeat' => 1,
'typecode' => 42
};
my $colname = $$data{colname};
push (@$colname, 'X', 'Y');
$$data{colname} = [ @$colname ];
writeFITSColumn(
file => $subSrcList,
extension => 'SRCLIST',
column => $data
);
info("Generated X and Y coordinate columns in $subSrcList");
# ---- completes task filter_list.pro
# ---- starting mask.pro
#
# ++++ Source counts VS mask radius lookup table
#
# Copy lookup table to intermediate file
# to have processing record of values used
# Should this be a CCF of some sort?
# Define lookup table intermediate
my $lookupFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Counts radius lookup table'
, format => 'ASCII'
) or return exception();
my $lookupInput = $ENV{PCMS_ROOT} . '/pipeline/etc/LOOKUP.FITS';
copyFile(
source => $lookupInput
, destination => $lookupFile
) or return exception();
info("Copied mask radius lookup table in place");
# Read FITS lookup table
my $lookupData = readFITSTable(
file => $lookupFile
, extension => 'LOOKUP'
, colname => [qw(CTS_LOW CTS_HIGH RADIUS OFFA_LOW OFFA_HIGH)]
) or return exception();
# Get number of rows from FITS table
my $num_lookup = numberFITSRows(
file => $lookupFile
, extension => 'LOOKUP'
);
# Reform lookup table from hash of arrays to array of hashes
my @lookup;
for ($i = 0; $i < $num_lookup; $i++) {
push(@lookup, {
'cts_low' => $lookupData->{CTS_LOW}[$i],
'cts_high' => $lookupData->{CTS_HIGH}[$i],
'radius' => $lookupData->{RADIUS}[$i],
'offa_low' => $lookupData->{OFFA_LOW}[$i],
'offa_high' => $lookupData->{OFFA_HIGH}[$i]
});
}
info("Read lookup table $lookupFile...");
# Get the OFFAX maximum value from $lookupFile = @lookup
my @offax_;
for $i ( @lookup ) {
push(@offax_ , $i->{offa_high} );
}
my $offax_lookupMax = max(@offax_);
# Create blank column arrays for mask radius
my (@allR_raw, @allR_ima);
# Loop over number of sources
for ($i = 0; $i < $num_sources; $i++) {
my $cts = $$mlCtsCol[$i]; # Counts in camera for this source
my $offax = $$mlOffaxCol[$i]; # offaxis angle for this source
# Limit OFFAX to maximum OFFAX in $lookupFile = @lookup
if ( $offax > $offax_lookupMax ){
my $new_offax = $offax_lookupMax - 0.001 ;
info("Source $i: OFFAX value $offax > max value $offax_lookupMax. Clamped OFFAX = $new_offax");
$offax = $new_offax ;
}
my @result = grep { $_->{cts_low} <= $cts && $_->{cts_high} > $cts &&
$_->{offa_low} <= $offax && $_->{offa_high} > $offax } @lookup;
my $rad = ($result[0]->{'radius'});
exception("Radius undefined for CTS $cts OFFAX $offax") unless defined($rad);
info("Source $i: CTS $cts OFFAX $offax MASK $rad");
# Write lookup radius to output column
$allR_raw[$i] = $rad*20;
$allR_ima[$i] = $rad/4;
}
# Get WCS reference keywords from header
my $refx1 = readFITSKeyword(
file => $epicImage
, extension => 'PRIMARY'
, keyword => 'REFXCRPX'
) or return exception();
my $refy1 = readFITSKeyword(
file => $epicImage
, extension => 'PRIMARY'
, keyword => 'REFYCRPX'
) or return exception();
my $crpix1 = readFITSKeyword(
file => $epicImage
, extension => 'PRIMARY'
, keyword => 'CRPIX1'
) or return exception();
my $crpix2 = readFITSKeyword(
file => $epicImage
, extension => 'PRIMARY'
, keyword => 'CRPIX2'
) or return exception();
# Loop over number of sources in source list
my (@allX_ima, @allY_ima);
# Convert from raw to image coordinates
for ($i = 0; $i < $num_sources; $i++) {
$allX_ima[$i] = (($allX_raw[$i] - $refx1)/80.0) + $crpix1;
$allY_ima[$i] = (($allY_raw[$i] - $refy1)/80.0) + $crpix2;
}
# Add new columns to source list
# rad_src_raw, xima, yima, rad_src_ima
$$data{RAD_RAW} = [@allR_raw];
$$data{RAD_IMA} = [@allR_ima];
$$data{X_IMA} = [@allX_ima];
$$data{Y_IMA} = [@allY_ima];
$$data{-column}{RAD_RAW} = { 'width' => 4,
'ttype' => 'RAD_RAW',
'repeat' => 1,
'typecode' => 42
};
$$data{-column}{RAD_IMA} = { 'width' => 4,
'ttype' => 'RAD_IMA',
'repeat' => 1,
'typecode' => 42
};
$$data{-column}{X_IMA} = { 'width' => 4,
'ttype' => 'X_IMA',
'repeat' => 1,
'typecode' => 42
};
$$data{-column}{Y_IMA} = { 'width' => 4,
'ttype' => 'Y_IMA',
'repeat' => 1,
'typecode' => 42
};
my $colname = $$data{colname};
pop @$colname;
pop @$colname;
push (@$colname, 'RAD_RAW', 'RAD_IMA', 'X_IMA', 'Y_IMA');
$$data{colname} = [ @$colname ];
writeFITSColumn(
file => $subSrcList,
extension => 'SRCLIST',
column => $data
);
info("Added RAD_RAW, RAD_IMA, X_IMA, Y_IMA columns to $subSrcList");
# FITS region file for emask
my $extrRegion = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'extraction region file'
, format => 'FITS'
) or return exception();
my $extrTemp = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'extraction temp file'
, format => 'FITS'
) or return exception();
doCommand(
'fcopy'
, infile => $subSrcList."[1][col X=X_IMA;Y=Y_IMA]"
, outfile => $extrTemp
) or return exception();
doCommand(
'fmerge'
, infiles => $extrTemp
, outfile => $extrRegion
, columns => "X,Y"
, mextname => 'REGION'
) or return exception();
# create shape column
my @SHAPEcol = ("CIRCLE") x $num_sources;
# write column to file
my $data2;
$$data2{SHAPE} = [@SHAPEcol];
$$data2{R} = [@allR_ima];
$$data2{-column}{SHAPE} = { 'width' => 10,
'colnum' => 1,
'ttype' => 'SHAPE',
'repeat' => 10,
'typecode' => 16
};
$$data2{-column}{R} = { 'width' => 4,
'ttype' => 'R',
'repeat' => 1,
'typecode' => 42
};
my $colname2;
push (@$colname2, 'R','SHAPE');
$$data2{colname} = [ @$colname2 ];
writeFITSColumn(
file => $extrRegion,
extension => 'REGION',
column => $data2
);
info("Generated region file $extrRegion");
# Generate new extraction detmask from region file
# using emask task and exposure image
my $extrDetMask = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'extraction detmask'
, format => 'FITS'
) or return exception();
my $expImage = findFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'EPIC exposure map'
, format => 'FITS'
) or return exception();
doCommand(
'emask'
, expimageset => $expImage
, detmaskset => $extrDetMask
, withregionset => 'yes'
, regionset => $extrRegion
) or return exception();
info("Generated extraction detmask $extrDetMask");
# check for RGA/OOT mask and combine with extraction detmask
my $pileupMask = findFile(
class => 'intermediate'
, instrument => thisInstrument
, content => 'Pileup mask'
, band => 8
, format => 'FITS'
);
if (fileExists( file => $pileupMask ) ) {
info("Combine extraction detmask with RGA/OOT mask for source extraction");
doCommand(
'farith'
, infil1 => "${extrDetMask}[1]"
, infil2 => $pileupMask
, outfil => $extrDetMask
, clobber => 'yes'
, copyprime => 'yes'
, ops => 'MUL'
) or return exception();
}
# ---- end mask.pro
# generate X,Y for SSP source list (use corrected RA, Dec)
my (@X_raw, @Y_raw);
my (@X_ima, @Y_ima);
my (@sCCD);
for ($j = 0; $j < $num_ssp; $j++) {
my $sc = $j+1;
my ($xcoord, $ycoord, $sccd) = coordConv($epicImage, $$srcCorrRaCol[$j], $$srcCorrDecCol[$j]);
$X_raw[$j] = int($xcoord);
$Y_raw[$j] = int($ycoord);
$X_ima[$j] = (($X_raw[$j] - $refx1)/80.0) + $crpix1;
$Y_ima[$j] = (($Y_raw[$j] - $refy1)/80.0) + $crpix2;
$sCCD[$j] = int($sccd);
info("DEBUG: Source $sc of $num_ssp SSPs: X $X_raw[$j] Y: $Y_raw[$j] CCD: $sCCD[$j]");
}
# generate R for SSP source list
my (@R_raw, @R_ima);
my @R_arcsec; # R in arcsec
my @R1_arcsec; # R1 in arcsec
for ($j = 0; $j < $num_ssp; $j++ ) {
my $cts = $$srcCtsCol[$j]; # counts in camera for this source
my $offax = $srcOffaxCol[$j]; # offaxis angle for this source
# Limit OFFAX to maximum OFFAX in $lookupFile = @lookup
if ( $offax > $offax_lookupMax ){
my $new_offax = $offax_lookupMax - 0.001 ;
info("Source $i: OFFAX value $offax > max value $offax_lookupMax. Clamped OFFAX = $new_offax");
$offax = $new_offax ;
}
my @result = grep { $_->{cts_low} < $cts && $_->{cts_high} > $cts &&
$_->{offa_low} < $offax && $_->{offa_high} > $offax } @lookup;
my $rad = ($result[0]->{'radius'});
exception("Radius undefined for CTS $cts OFFAX $offax") unless defined($rad);
# write raw, image radius to output column
$R_raw[$j] = $rad * 20.0;
$R_ima[$j] = $rad / 4.0;
}
# ---- return to spec_noregion.pro
# Loop for each bright source (in SSP list), select circular region
# with >= 90% good area for background extraction (count pixels
# in region in new detector mask).
# Start first background loop
info("==== Creating background regions (step 1)");
my $MAX_M = 150; # Maximum iteration number
# Background position and radius arrays
# declare these here because they will be used later
## my (@XB_raw, @YB_raw, @RB_raw);
my (@BKG_SHAPE, @XB_raw, @YB_raw, @RB_raw, @RBi_raw, @RBo_raw);
# Switches for PN background Yes or Not in case no Bkg is possible for some source.
my @BKG;
my $good_area_limit = 0.70;
my $illum = 0.0;
my $max_illum = 0.0;
# Get observation submode
my $submode = readFITSKeyword(
file => $forSpecFilteredEvList
, extension => 'EVENTS'
, keyword => 'SUBMODE'
) or return exception();
info("DEBUG: Observation submode is $submode");
#
# No iterations since 'ebkgreg' task for background region estimation
#
# my $smallWindowMode = 0;
# if (thisInstrument =~ /emos/i && $submode =~/PrimePartialW/i) {
# $smallWindowMode = 1;
# info("Small Window Mode $cam:$submode detected, changing iteration parameters");
# }
# Some variables that will be used for forked parallelization
my $nforks = $ENV{PCMS_NFORKS};
my $jobrangesref;
my $usedforks;
my @children; # child process IDs
my $badexitstatus;
($jobrangesref,$usedforks) = jobRanges( $num_ssp, $nforks );
Pcmsdb::set_InactiveDestroy(1); # So that children won't kill DB connections
foreach my $job (0..$usedforks-1) {
my $pid = fork();
if ($pid) {
# this is a parent, just keep track of child IDs
push @children, $pid;
} elsif ($pid == 0) {
eval {
# this is a child. Here is where the work happens.
my $range1 = ${$jobrangesref}[$job];
my $range2 = ${$jobrangesref}[$job+1]-1;
info("********** Start of forked child process $$ for sources $range1 -- $range2");
# These are only used internally by the forked process
my @RB_ima = @allR_ima; # Set bg radius = source radius initially
my (@XB_ima, @YB_ima);
# Loop over number of bright sources
for ($j = $range1; $j <= $range2; $j++) {
info("--------------------------------------------------------------");
my $id = $$srcIdCol[$j];
my $sc = $j+1;
# Check if source is faint in this camera
my $camcts = $$srcCtsCol[$j];
info("**** Source $id ($sc of $num_ssp): In-camera counts $camcts, CCD: $sCCD[$j]");
info("**** Source $id ($sc of $num_ssp): Step 1 background");
####------------------------------------------------------------------------------------------------
### # Background selection method for MOS's
### if ((thisInstrument eq 'emos1' || thisInstrument eq 'emos2')) {
###
### # Set parameters for iteration
### # Increase iterator only for SmallWindow mode sources in CCD 1
###
### my $m = ($smallWindowMode && $sCCD[$j] == 1) ? 45 : 0; # set iterator value
###
### # Start optimisation loop over offset parameter m
###
### while ($m < $MAX_M) {
###
### $m = $m + 1;
### info("[1] Source ID: $id, iteration $m: Running select_bkg(".$X_ima[$j].",".$Y_ima[$j].",".$R_ima[$j].")");
###
### my $result = &select_bkg($X_ima[$j], $Y_ima[$j], $R_ima[$j], $m, $extrDetMask, $epicImage, $id);
###
### # Get iteration illumination value
###
### $illum = $$result->{illum};
### my $max_index = $$result->{index};
###
### # Illumination exceeds good illumination limit
###
### if ($illum > $good_area_limit) {
###
### info("Source $id: CONVERGED @ radial step: $m azimuthal step: $max_index illumination: $illum");
###
### # Get optimised values of XB, YB, RB in image coordinates
###
### $XB_ima[$j] = $$result->{xback};
### $YB_ima[$j] = $$result->{yback};
### $RB_ima[$j] = $$result->{radback};
###
### info("Source: $id: background solution XB_ima, YB_ima, RB_ima: ".$XB_ima[$j] . "," . $YB_ima[$j] . "," . $RB_ima[$j]);
###
### last; # halt the while... loop
###
### # Illumination exceeds best value so far
###
### } elsif ($illum > $max_illum) {
###
### $XB_ima[$j] = $$result->{xback};
### $YB_ima[$j] = $$result->{yback};
### $RB_ima[$j] = $$result->{radback};
###
### $max_illum = $illum; # update current best illumination
### }
###
### }
###
### # Get raw values of xb, yb, rb
###
### $XB_raw[$j] = ($XB_ima[$j] - $crpix1) * 80.0 + $refx1;
### $YB_raw[$j] = ($YB_ima[$j] - $crpix2) * 80.0 + $refy1;
###
### $RB_raw[$j] = $RB_ima[$j] * 80.0;
###
### # Reset max illumination
### $max_illum = 0;
###
####------------------------------------------------------------------------------------------------
### # Background selection method for PN
### } elsif ( thisInstrument eq 'epn' ) {
#
# Same algorithm to find the best Bkg extraction region for PN and MOSs
#
# The algorith searchs for a circular background region in the same CCD where the source is located, except for the source in the central CCD
# of a MOS observation in SmallWindow mode (PrimePartialW2/3). In that case the background is estimated from an annulus (inner radius = 5.5 arcmin,
# outer radius = 11 arcmin) centered in the center of the image. Thus the background is estimated from the peripheral CCDs and the central CCD is
# completely excluded.
#
# For EPIC-pn sources the algorithm avoids the same RAWY column of the source in order to exclude out-of-time events from the background estimation.
#
# The background region always has a radius larger than 3 pixels, otherwise no background is calculated.
#
$R_arcsec[$j] = $R_ima[$j] * 4;
info("[1] Source ID: $id, Running select_bkg_ebkgreg(".$X_raw[$j].",".$Y_raw[$j].",".$R_arcsec[$j].")");
# $R_ima[$j] = $rad / 4.0 . LOOKUP -> $rad (arcsec) ==> ebkgreg( r = R_ima * 4 )
my $result_ebkgreg = &select_bkg_ebkgreg($X_raw[$j], $Y_raw[$j], $R_arcsec[$j], $epicImage, $id);
$BKG_SHAPE[$j] = $$result_ebkgreg->{shape};
$XB_raw[$j] = $$result_ebkgreg->{xback};
$YB_raw[$j] = $$result_ebkgreg->{yback};
$RBi_raw[$j] = $$result_ebkgreg->{radbackinner};
$RBo_raw[$j] = $$result_ebkgreg->{radbackouter};
$RB_raw[$j] = $RBi_raw[$j];
info("DEBUG 1st - ebkgreg - Source ID: $id = (". $BKG_SHAPE[$j] . "," .$XB_raw[$j].",".$YB_raw[$j].",".$RBi_raw[$j].",".$RBo_raw[$j].")");
if ( ! ( defined($XB_raw[$j]) && defined($YB_raw[$j]) && defined($RB_raw[$j]) ) ){ # <<<===========================================================================
### $PN_BKG[$j] = 'NONE';
### info("DEBUG 1st - ebkgreg - Source ID: $id = No background is possible. PN_BKG = $PN_BKG[$j]");
$BKG[$j] = 'NONE';
info("DEBUG 1st - ebkgreg - Source ID: $id = No background is possible. BKG = $BKG[$j]");
### } else { $PN_BKG[$j] = '';}
} else { $BKG[$j] = '';}
### }
#------------------------------------------------------------------------------------------------
} # end loop over sources
}; # end eval
# This is the end of the forked child process. Write the
# seqdb to a temporary file which we will use to
# communicate information back to the parent process
# (e.g., the names of new files that were created). If
# an exception was detected inside the eval {} block,
# we also record it to the seqdb file.
my $exception;
if ($@) {
$exception = $@;
}
# The ?B_raw variables are needed in the next section, so
# we will attach it to the child seqdb file, and then load
# them back in the parent.
writeChildSeqdb(exception => $exception,
extra => { XB_raw_ref => \@XB_raw,
YB_raw_ref => \@YB_raw,
RB_raw_ref => \@RB_raw,
BKG_ref => \@BKG, # Before: PN_BKG_ref => \@PN_BKG,
BKG_SHAPE_ref => \@BKG_SHAPE } # <<<==============================================================================================
);
# Return status of child process reflects occurence of an exception
if ($exception) {
exit -1;
} else {
exit 0;
}
} else {
# Couldn't fork
return exception("Couldn't fork: $!\n");
}
} # end loop over jobs
# Now wait for all the child jobs to finish.
# Check to see if any gave a bad exit status before returning.
$badexitstatus = waitChild( \@children );
# Necessary so that DB connections may be destroyed again by the parent
Pcmsdb::set_InactiveDestroy(0);
# Merge sequence information from the children
my $extraarrayref = mergeChildSeqdb(\@children);
if ($badexitstatus) {
# mergeChildSeqdb should have copied any exception in a child
# process into the parent Seqdb which will be picked up in the
# following exception() call.
return exception();
}
# If the children were OK, extract extra information
foreach my $job(0..$usedforks-1) {
my $range1 = ${$jobrangesref}[$job];
my $range2 = ${$jobrangesref}[$job+1]-1;
my $extraref = ${$extraarrayref}[$job];
# The [X/Y/R]B_raw arrays from each child were stored in the seqdb file
@BKG_SHAPE[$range1..$range2] = @{${$extraref}{BKG_SHAPE_ref}}[$range1..$range2];
@XB_raw[$range1..$range2] = @{${$extraref}{XB_raw_ref}}[$range1..$range2];
@YB_raw[$range1..$range2] = @{${$extraref}{YB_raw_ref}}[$range1..$range2];
@RB_raw[$range1..$range2] = @{${$extraref}{RB_raw_ref}}[$range1..$range2];
### @PN_BKG[$range1..$range2] = @{${$extraref}{PN_BKG_ref}}[$range1..$range2]; # <<<=======================================================================
@BKG[$range1..$range2] = @{${$extraref}{BKG_ref}}[$range1..$range2];
}
info("==== Step 1 iteration complete!");
# STEP 2: Optimise selection regions with eregionanalyse
info("==== Step 2 eregionanalyse optimisation");
# The following arrays will be used in later steps
my (@X1_raw, @Y1_raw, @R1_raw);
my (@X1_ima, @Y1_ima, @R1_ima);
my @source_radius = ();
($jobrangesref,$usedforks) = jobRanges( $num_ssp, $nforks );
Pcmsdb::set_InactiveDestroy(1); # So that children won't kill DB connections
undef @children;
foreach my $job (0..$usedforks-1) {
my $pid = fork();
if ($pid) {
# this is a parent, just keep track of child IDs
push @children, $pid;
} elsif ($pid == 0) {
eval {
# this is a child. Here is where the work happens.
my $range1 = ${$jobrangesref}[$job];
my $range2 = ${$jobrangesref}[$job+1]-1;
info("********** Start of forked child process $$ for sources $range1 -- $range2");
for ( $j = $range1; $j <= $range2; $j++ ) {
info("------------------------------------------------------------");
my $id = $$srcIdCol[$j];
my $sc = $j+1;
info("++++ Source ID $id ($sc of $num_ssp): Starting eregionanalyse");
########################### <<============================================================================================
#
# No eregionanalyse for sources for which Bkg region was not possible
#
### if ( ( thisInstrument eq 'epn' ) && ( $PN_BKG[$j] eq 'NONE' ) ){
### info("Background in PN was not possible for Source: $id, skipping eregionanalyse.");
### next;
### }
my $instrument = thisInstrument;
if ( $BKG[$j] eq 'NONE' ){
info("Background in $instrument was not possible for Source: $id, skipping eregionanalyse.");
next;
}
########################### <<============================================================================================
# This loop ('eregionanalyse') only runs for Bkg region shape = CIRCLE. For ANNULUS it is skipped:
if ( $BKG_SHAPE[$j] eq 'CIRCLE' ){
my $cmd = 'eregionanalyse centroid=no imageset='.$epicImage.' srcexp="(X,Y) IN circle(' .
"$X_raw[$j], $Y_raw[$j], $R_raw[$j])" . '" backexp="(X,Y) IN circle(' .
"$XB_raw[$j], $YB_raw[$j], $RB_raw[$j])" . '"';
info("CMD: $cmd");
my $optoutput = `$cmd`;
info("OPTOUTPUT: $optoutput");
# Get optimised source X', Y', R' from eregionanalyse output
my ($counts) = ($optoutput =~ /counts\sin\ssource\sregion:\s+(\d+)/);
info("Found $counts counts in region");
if ($counts == 0) {
# Ignore optimisation in case of zero counts in region
info("Ignoring optimisation...");
$X1_raw[$j] = $X_raw[$j];
$Y1_raw[$j] = $Y_raw[$j];
$R1_raw[$j] = $R_raw[$j];
} else {
($X1_raw[$j]) = ($optoutput =~ /xcentroid:\s+([\d.]+)/);
($Y1_raw[$j]) = ($optoutput =~ /ycentroid:\s+([\d.]+)/);
($R1_raw[$j]) = ($optoutput =~ /optradius:\s+\d+\s+arcsecs\s+(\d+)/);
}
if ($R1_raw[$j] < 240.0) {$R1_raw[$j] = 240.0}; # clamp optimised radius to minimum 240
if ($R1_raw[$j] > 800.0) {$R1_raw[$j] = 800.0}; # clamp optimised radius to maximum 800
} elsif ( $BKG_SHAPE[$j] eq 'ANNULUS' ){
info("Background region for $instrument is $BKG_SHAPE[$j] for Source: $id, skipping eregionanalyse.");
# Optimized Src position = Previous source position
info("Ignoring optimisation with eregionanalyse...");
$X1_raw[$j] = $X_raw[$j];
$Y1_raw[$j] = $Y_raw[$j];
$R1_raw[$j] = $R_raw[$j];
}
# convert optimised source X', Y', R' to image coordinates
$X1_ima[$j] = (($X1_raw[$j] - $refx1)/80.0) + $crpix1;
$Y1_ima[$j] = (($Y1_raw[$j] - $refy1)/80.0) + $crpix2;
$R1_ima[$j] = $R1_raw[$j] / 80.0;
# save copy of extraction radius indexed by source ID
$source_radius[$id] = $R1_raw[$j];
} # end loop over sources
}; # end eval
# end of the child process
my $exception;
if ($@) {
$exception = $@;
}
writeChildSeqdb(exception => $exception,
extra => { X1_raw_ref => \@X1_raw,
Y1_raw_ref => \@Y1_raw,
R1_raw_ref => \@R1_raw,
X1_ima_ref => \@X1_ima,
Y1_ima_ref => \@Y1_ima,
R1_ima_ref => \@R1_ima,
source_radius_ref => \@source_radius }
);
# Return status of child process reflects occurence of an exception
if ($exception) {
exit -1;
} else {
exit 0;
}
} else {
# Couldn't fork
return exception("Couldn't fork: $!\n");
}
} # end loop over jobs
# Now wait for all the child jobs to finish.
# Check to see if any gave a bad exit status before returning.
$badexitstatus = waitChild( \@children );
# Necessary so that DB connections may be destroyed again by the parent
Pcmsdb::set_InactiveDestroy(0);
# Merge sequence information from the children
my $extraarrayref = mergeChildSeqdb(\@children);
if ($badexitstatus) {
# mergeChildSeqdb should have copied any exception in a child
# process into the parent Seqdb which will be picked up in the
# following exception() call.
return exception();
}
# If the children were OK, extract extra information
foreach my $job(0..$usedforks-1) {
my $range1 = ${$jobrangesref}[$job];
my $range2 = ${$jobrangesref}[$job+1]-1;
my $extraref = ${$extraarrayref}[$job];
# The [X/Y/R]1_[raw/ima], source_radius arrays from each child were stored in the seqdb file
@X1_raw[$range1..$range2] = @{${$extraref}{X1_raw_ref}}[$range1..$range2];
@Y1_raw[$range1..$range2] = @{${$extraref}{Y1_raw_ref}}[$range1..$range2];
@R1_raw[$range1..$range2] = @{${$extraref}{R1_raw_ref}}[$range1..$range2];
@X1_ima[$range1..$range2] = @{${$extraref}{X1_ima_ref}}[$range1..$range2];
@Y1_ima[$range1..$range2] = @{${$extraref}{Y1_ima_ref}}[$range1..$range2];
@R1_ima[$range1..$range2] = @{${$extraref}{R1_ima_ref}}[$range1..$range2];
foreach my $j ($range1..$range2) {
my $id = $$srcIdCol[$j];
$source_radius[$id] = ${${$extraref}{source_radius_ref}}[$id];
}
}
# Diagnostic loop
info("*** DEBUG DUMP");
for ($j=0; $j < $num_ssp; $j++) {
info("Source $j: raw (".$X1_raw[$j].",".$Y1_raw[$j].",".$R1_raw[$j].") ima (".$X1_ima[$j].",".$Y1_ima[$j].",".$R1_ima[$j].")");
}
info("*** DEBUG DUMP");
# Start second background optimisation loop
# Loop over number of sources
info("==== Step 3: Creating background");
my $MAX_M = 150; # Maximum iteration number
# These will be used in a later section
my (@XB1_raw, @YB1_raw, @RB1_raw, @RB1i_raw, @RB1o_raw);
($jobrangesref,$usedforks) = jobRanges( $num_ssp, $nforks );
Pcmsdb::set_InactiveDestroy(1); # So that children won't kill DB connections
undef @children;
foreach my $job (0..$usedforks-1) {
my $pid = fork();
if ($pid) {
# this is a parent, just keep track of child IDs
push @children, $pid;
} elsif ($pid == 0) {
eval {
# this is a child. Here is where the work happens.
my $range1 = ${$jobrangesref}[$job];
my $range2 = ${$jobrangesref}[$job+1]-1;
info("********** Start of forked child process $$ for sources $range1 -- $range2");
# These are only used locally
my (@XB1_ima, @YB1_ima, @RB1_ima);
for ( $j = $range1; $j <= $range2; $j++ ) {
my $id = $$srcIdCol[$j];
my $sc = $j+1;
info("**** Source ID $id ($sc of $num_ssp): Step 3 background");
########################### <<============================================================================================
#
# No 2nd ebkgreg for sources for which Bkg region was not possible
#
### if ( ( thisInstrument eq 'epn' ) && ( $PN_BKG[$j] eq 'NONE' ) ){
### info("Background in PN was not possible for Source: $id, skipping 2nd ebkgreg run.");
### next;
### }
my $instrument = thisInstrument;
if ( $BKG[$j] eq 'NONE' ){
info("Background in $instrument was not possible for Source: $id, skipping 2nd ebkgreg run.");
next;
}
########################### <<============================================================================================
####------------------------------------------------------------------------------------------------
### # Background selection method for MOS's
### if ((thisInstrument eq 'emos1' || thisInstrument eq 'emos2')) {
###
### if ( $X1_ima[$j] > 0.0 && $Y1_ima[$j] > 0.0 ) { # Check for non-negative X', Y'
###
### # Set parameters for iteration
### # Increase iterator only for SmallWindow mode sources in CCD 1
###
### my $m = ($smallWindowMode && $sCCD[$j] == 1) ? 45 : 0; # set iterator value
###
### while ($m < $MAX_M) {
###
### $m = $m + 1;
### info("[3] Source ID: $id, iteration $m: Running select_bkg(".$X1_ima[$j].",".$Y1_ima[$j].",".$R1_ima[$j].")");
###
### my $result = &select_bkg($X1_ima[$j], $Y1_ima[$j], $R1_ima[$j], $m, $extrDetMask, $epicImage, $id);
###
### # get maximum illumination value
###
### $illum = $$result->{illum};
### my $max_index = $$result->{index};
###
### # Illumination exceeds min 'good' value
###
### if ($illum > $good_area_limit) {
###
### info("Source $id: CONVERGED @ radial step: $m azimuthal step: $max_index illumination: $illum");
###
### # Get optimised values of XB, YB, RB in image coordinates
###
### $XB1_ima[$j] = $$result->{xback};
### $YB1_ima[$j] = $$result->{yback};
### $RB1_ima[$j] = $$result->{radback};
###
### info("Source: $id: background solution XB1_ima, YB1_ima, RB1_ima: ".$XB1_ima[$j] . "," . $YB1_ima[$j] . "," . $RB1_ima[$j]);
###
### last; # Halt the while... loop
###
### # Illumination exceeds best value so far
###
### } elsif ($illum > $max_illum) {
###
### $XB1_ima[$j] = $$result->{xback};
### $YB1_ima[$j] = $$result->{yback};
### $RB1_ima[$j] = $$result->{radback};
###
### $max_illum = $illum; # update current best illumination
### }
###
### }
###
### }
###
### # Get raw values for xb', yb', rb'
###
### $XB1_raw[$j] = ($XB1_ima[$j] - $crpix1) * 80.0 + $refx1;
### $YB1_raw[$j] = ($YB1_ima[$j] - $crpix2) * 80.0 + $refy1;
###
### $RB1_raw[$j] = $RB1_ima[$j] * 80.0;
###
### # Reset max illumination
###
### $max_illum = 0;
###
####------------------------------------------------------------------------------------------------
### # Background selection method for PN
### } elsif ( thisInstrument eq 'epn' ) {
$R1_arcsec[$j] = $R1_ima[$j] * 4;
info("[3] Source ID: $id, Running select_bkg_ebkgreg(".$X1_raw[$j].",".$Y1_raw[$j].",".$R1_arcsec[$j].")");
# $R_ima[$j] = $rad / 4.0 . LOOKUP -> $rad (arcsec) ==> ebkgreg( r = R_ima * 4 )
my $result_ebkgreg = &select_bkg_ebkgreg($X1_raw[$j], $Y1_raw[$j], $R1_arcsec[$j], $epicImage, $id);
$BKG_SHAPE[$j] = $$result_ebkgreg->{shape};
$XB1_raw[$j] = $$result_ebkgreg->{xback};
$YB1_raw[$j] = $$result_ebkgreg->{yback};
$RB1i_raw[$j] = $$result_ebkgreg->{radbackinner};
$RB1o_raw[$j] = $$result_ebkgreg->{radbackouter};
$RB1_raw[$j] = $RB1i_raw[$j];
info("DEBUG 2nd - ebkgreg - Source ID: $id = (". $BKG_SHAPE[$j] . "," .$XB1_raw[$j].",".$YB1_raw[$j].",".$RB1i_raw[$j].",".$RB1o_raw[$j].")");
### }
#------------------------------------------------------------------------------------------------
} # end loop over sources
}; # end eval
# end of the child process
my $exception;
if ($@) {
$exception = $@;
}
writeChildSeqdb(exception => $exception,
extra => { XB1_raw_ref => \@XB1_raw,
YB1_raw_ref => \@YB1_raw,
RB1i_raw_ref => \@RB1i_raw,
RB1o_raw_ref => \@RB1o_raw,
BKG_SHAPE_ref => \@BKG_SHAPE }
);
# Return status of child process reflects occurence of an exception
if ($exception) {
exit -1;
} else {
exit 0;
}
} else {
# Couldn't fork
return exception("Couldn't fork: $!\n");
}
} # end loop over jobs
# Now wait for all the child jobs to finish.
# Check to see if any gave a bad exit status before returning.
$badexitstatus = waitChild( \@children );
# Necessary so that DB connections may be destroyed again by the parent
Pcmsdb::set_InactiveDestroy(0);
# Merge sequence information from the children
my $extraarrayref = mergeChildSeqdb(\@children);
if ($badexitstatus) {
# mergeChildSeqdb should have copied any exception in a child
# process into the parent Seqdb which will be picked up in the
# following exception() call.
return exception();
}
# If the children were OK, extract extra information
foreach my $job(0..$usedforks-1) {
my $range1 = ${$jobrangesref}[$job];
my $range2 = ${$jobrangesref}[$job+1]-1;
my $extraref = ${$extraarrayref}[$job];
# The [X/Y/R]1_[raw/ima], source_radius arrays from each child were stored in the seqdb file
@XB1_raw[$range1..$range2] = @{${$extraref}{XB1_raw_ref}}[$range1..$range2];
@YB1_raw[$range1..$range2] = @{${$extraref}{YB1_raw_ref}}[$range1..$range2];
@RB1i_raw[$range1..$range2] = @{${$extraref}{RB1i_raw_ref}}[$range1..$range2];
@RB1o_raw[$range1..$range2] = @{${$extraref}{RB1o_raw_ref}}[$range1..$range2];
@BKG_SHAPE[$range1..$range2] = @{${$extraref}{BKG_SHAPE_ref}}[$range1..$range2];
}
info("==== Step 3 iteration complete");
# Iterate over sources and remove contaminating sources from bkg
# Generate region files for light curve and spectral extraction
info("==== Step 4: Start final extraction region generation");
($jobrangesref,$usedforks) = jobRanges( $num_ssp, $nforks );
Pcmsdb::set_InactiveDestroy(1); # So that children won't kill DB connections
undef @children;
foreach my $job (0..$usedforks-1) {
my $pid = fork();
if ($pid) {
# this is a parent, just keep track of child IDs
push @children, $pid;
} elsif ($pid == 0) {
eval {
# this is a child. Here is where the work happens.
my $range1 = ${$jobrangesref}[$job];
my $range2 = ${$jobrangesref}[$job+1]-1;
info("********** Start of forked child process $$ for sources $range1 -- $range2");
my @regions;
for ( $j = $range1; $j <= $range2; $j++ ) {
my $id = $$srcIdCol[$j];
my $snum = $j+1;
info("* * * Writing region for source ID: $id ($snum of $num_ssp)");
########################### <<============================================================================================
#
# No Src and Bkg regions for sources for which Bkg region was not possible
#
### if ( ( thisInstrument eq 'epn' ) && ( $PN_BKG[$j] eq 'NONE' ) ){
### info("Background in PN was not possible for Source: $id, skipping this source for this exposure.");
### next;
### }
my $instrument = thisInstrument;
if ( $BKG[$j] eq 'NONE' ){
info("Background in $instrument was not possible for Source: $id, skipping this source for this exposure.");
next;
}
########################### <<============================================================================================
# define source extraction region
my $sx = $X1_raw[$j];
my $sy = $Y1_raw[$j];
my $sr = $R1_raw[$j];
###########################
info("DEBUG - Source id: $id Shape: $BKG_SHAPE[$j] Source extraction region sx, sy, sr = $sx, $sy, $sr");
###########################
# generate source region files
my $srcFITSRegionFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'Source extraction region file'
, format => 'FITS'
);
writeRegionRow2($srcFITSRegionFile, "CIRCLE", $sx, $sy, $sr);
my $srcASCRegionFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'src ds9 region'
, format => 'ASCII'
);
my $rtn = writeASCIIRegion( $srcASCRegionFile
, "circle"
, $sx
, $sy
, $sr
, 0 # R_outer = 0 for circle
, "green"
, "$id"
, 0 );
# define background extraction region
my $bkg_shape = $BKG_SHAPE[$j];
my $bx = $XB1_raw[$j];
my $by = $YB1_raw[$j];
### my $br = $RB1_raw[$j];
my $br_i = $RB1i_raw[$j];
my $br_o = $RB1o_raw[$j];
my $br = $br_i;
# generate background region file
my $bkgFITSRegionFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'Background extraction region file'
, format => 'FITS'
);
my $bkgASCRegionFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'bkg ds9 region'
, format => 'ASCII'
);
### writeRegionRow($bkgFITSRegionFile, "CIRCLE", $bx, $by, $br);
writeRegionRow2($bkgFITSRegionFile, $bkg_shape, $bx, $by, $br_i, $br_o);
my $rtn = writeASCIIRegion( $bkgASCRegionFile
, $bkg_shape
, $bx
, $by
, $br_i
, $br_o # R_outer = 0 for circle
, "white"
, "$id"
, 0 ); # append = 0
# loop over all sources, find contamination
my $k;
for ($k=0; $k < $num_sources; $k++) {
if ( $id == $$mlIdCol[$k] ) { next; }
my $cx = $allX_raw[$k];
my $cy = $allY_raw[$k];
my $cr = $allR_raw[$k];
my $sep = sqrt(($bx - $cx)**2 +
($by - $cy)**2);
# check angular separation, modify region if found
if ( $bkg_shape =~ /CIRCLE/i) {
if ($sep < ($br + $cr)) {
info("Contaminating source at $cx, $cy");
# write exclusion row to corresponding background region file
writeRegionRow2($bkgFITSRegionFile, "!CIRCLE", $cx, $cy, $cr);
# append exclusion row to ASCII background region file
my $rtn = writeASCIIRegion( $bkgASCRegionFile
, "circle"
, $cx
, $cy
, $cr
, 0 # R_outer = 0 for circle
, "red"
, "$$mlIdCol[$k]"
, 1 ); # append = 1
} } elsif ( $bkg_shape =~ /ANNULUS/i ){
if ( $sep < ($br_o + $cr) && $sep > ($br_i - $cr) ) {
info("Contaminating source at $cx, $cy");
# write exclusion row to corresponding background region file
writeRegionRow2($bkgFITSRegionFile, "!CIRCLE", $cx, $cy, $cr);
# append exclusion row to ASCII background region file
my $rtn = writeASCIIRegion( $bkgASCRegionFile
, "circle"
, $cx
, $cy
, $cr
, 0 # R_outer = 0 for circle
, "red"
, "$$mlIdCol[$k]"
, 1 ); # append = 1
} }
}
# push source, background extraction region filenames onto arraay of hashes
push (@regions, { 'id' => $id,
'source' => $srcFITSRegionFile,
'background' => $bkgFITSRegionFile
});
} # end loop over sources
}; # end eval
# end of the child process
my $exception;
if ($@) {
$exception = $@;
}
writeChildSeqdb(exception => $exception);
# Return status of child process reflects occurence of an exception
if ($exception) {
exit -1;
} else {
exit 0;
}
} else {
# Couldn't fork
return exception("Couldn't fork: $!\n");
}
} # end loop over jobs
# Now wait for all the child jobs to finish.
# Check to see if any gave a bad exit status before returning.
$badexitstatus = waitChild( \@children );
# Necessary so that DB connections may be destroyed again by the parent
Pcmsdb::set_InactiveDestroy(0);
# Merge sequence information from the children
my $extraarrayref = mergeChildSeqdb(\@children);
if ($badexitstatus) {
# mergeChildSeqdb should have copied any exception in a child
# process into the parent Seqdb which will be picked up in the
# following exception() call.
return exception();
}
info("==== Step 4: Final extraction region generation complete");
# --------------------------------------------------------------------------------
##
## END: 2006-03-15 , DJF
##
# Sort the EXPOSU extensions of event lists.
# This has to be done for both PN and MOS
my $tempEventList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Temporary event list'
);
copyFile(source => $forTimeFilteredEvList, destination => $tempEventList);
foreach my $i (@numCcds)
{
my $extname = sprintf('EXPOSU%02d', $i);
if (!hasFITSExtension(file => $tempEventList, extension => $extname))
{
info("File $tempEventList doesn't contain extension $extname - sort won't be performed.");
next;
}
doCommand(
'fsort'
, infile => $tempEventList."[$extname]"
, columns => 'TIME'
, method => 'insert'
) or return exception();
}
copyFile(source => $tempEventList, destination => $forTimeFilteredEvList);
# Now make the products:
# Record the fact this exposure is used for source specific products in the database
setExposureProperty( instrument => thisInstrument , exp_id => $exp_id , name => 'ssp' , value => 0+TRUE );
my $allSrcRegionFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'all src ds9 region'
, format => 'ASCII'
);
doCommand(
'slconv'
, srclisttab => "${filteredMlList}:SRCLIST"
, radiusexpression => 60
, radiusstyle => 'raw'
, withlabels => 'no'
, colour => 'cyan'
, outputstyle => 'ds9'
, outfilestyle => 'whole'
, outfile => $allSrcRegionFile
)
or return exception();
# Parallelize the next portion of the code. We will split up the
# work into "jobs" for each of N forked processed, where N should
# be related to the number of cores available. The logging system
# seems to be based on Log::Dispatch which claims to be
# fork/thread safe:
# 'This module will close and re-open the logfile after a fork.'
# While log messages may be confusing, they should not walk over
# eachother while actually writing a given message string to file.
# Figure out how to divide up the work amongst the child processes
my $nselectedSources = scalar @selectedSources;
($jobrangesref,$usedforks) = jobRanges( $nselectedSources, $nforks );
Pcmsdb::set_InactiveDestroy(1); # So that children won't kill DB connections
undef @children;
foreach my $job (0..$usedforks-1) {
my $pid = fork();
if ($pid) {
# this is a parent, just keep track of child IDs
push @children, $pid;
} elsif ($pid == 0) {
# We put the contents of the forked process into an
# eval. The reason for this is that doCommand can generate
# an exception which will exit the module
# immediately. This way we can catch the exceptions, and
# at the end of the forked portion (after the eval) we can
# first add the error to the Seqdb and write it to a file
# for this PID. The parent can then check the file for the
# exception. We also exit the child with nonzero status.
eval {
# this is a child. Here is where the work happens.
my $range1 = ${$jobrangesref}[$job];
my $range2 = ${$jobrangesref}[$job+1]-1;
info("********** Start of forked child process $$ for sources $range1 -- $range2");
foreach my $i (@selectedSources[$range1..$range2])
{
#if ($job == 0) {
# info("Generate a test exception using a false doCommand");
# doCommand(
# 'borkborkbork'
# ) or return exception();
#}
my @kwdNameList = qw(SRCNUM SRC_RA SRC_DEC SRC_CRA SRC_CDEC); # get FLAG here? GD
my @kwdTypeList = qw(integer real real real real);
my @kwdRealValsList = ($$srcRaCol[$i], $$srcDecCol[$i]
, $$srcCorrRaCol[$i], $$srcCorrDecCol[$i]);
my @kwdWithCommentList = qw(y y y y y);
my @kwdCommentList = ("EPIC source list number"
, "RA of source", "DEC of source", "Corrected RA of source"
, "Corrected DEC of source");
#### units not yet available as a parameter to addAttribute?
my $id=$$srcIdCol[$i];
my $src_ra=$$srcRaCol[$i];
my $src_dec=$$srcDecCol[$i];
my $row_in_filteredML=$mlId_to_rowNum{ $mlSrcIdColInSSP->[$i] };
########################### <<============================================================================================
#
# No Products ( LC and Spectra ) for sources for which Bkg region was not possible
#
### if ( ( thisInstrument eq 'epn' ) && ( $PN_BKG[$j] eq 'NONE' ) ){
### info("Background in PN was not possible for Source: $id, skipping product generation.");
### next;
### }
my $instrument = thisInstrument;
if ( $BKG[$i] eq 'NONE' ){
info("Background in $instrument was not possible for Source: $id, skipping product generation.");
next;
}
########################### <<============================================================================================
my $srcRegionSet = findFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Source extraction region file'
, format => 'FITS'
);
my $bkgRegionSet = findFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Background extraction region file'
, format => 'FITS'
);
# Check whether source lies within the detector window with Richard's new
# task 'echeckregion' in the especget package; should replace the subsequent
# check on the bins of the light curve (see below) but keep that one in as
# a safe guard
my $echeckregion_output = newFile(class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'echeckregion output'
, format => 'ASCII'
);
my $echeckregion_status = ModuleUtil::isSourceOnChip(
"region($srcRegionSet,X,Y)"
, $forTimeFilteredEvList
, $echeckregion_output
);
if (!$echeckregion_status)
{
info("Source $id lies outside detector window, skipping this source for this exposure.");
next;
} else {
info("Source $id seems to be within detector window. Continue making products.");
}
# Make and process the light curve for this source:
# x1 BINNING --------------------------------------
my $srcRateSet = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Uncorrected source time series'
);
my $bgdRateSet = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Uncorrected background time series'
);
doCommand(
'evselect'
, table => "${forTimeFilteredEvList}:EVENTS"
, writedss => 'yes'
, updateexposure => 'yes'
, expression => "region($srcRegionSet,X,Y)"
, energycolumn => 'PI'
, withrateset => 'yes'
, rateset => $srcRateSet
, timebinsize => $$timeBinSizeCol[$i]
, maketimecolumn => 'yes'
, makeratecolumn => 'yes'
)
or return exception();
doCommand(
'evselect'
, table => "${forTimeFilteredEvList}:EVENTS"
, writedss => 'yes'
, updateexposure => 'yes'
, expression => "region($bkgRegionSet,X,Y)"
, energycolumn => 'PI'
, withrateset => 'yes'
, rateset => $bgdRateSet
, timebinsize => $$timeBinSizeCol[$i]
, maketimecolumn => 'yes'
, makeratecolumn => 'yes'
)
or return exception();
# Check whether the lightcurve has any bins, otherwise skip any further
# steps. This catches cases where, eg. full-frame imaging overlaps window
# mode data (so sources exist in the list where there are no photons
# in the windowed exposure)
my @rate = readFITSColumn(file => $srcRateSet
, extension => 'RATE'
, column => 'RATE'
);
unless(grep { $_ } @rate)
{
info("Lightcurve contains no useful bins, skipping this source for this exposure");
next;
}
my $bkgSubtractedRateSet = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Exposure-corrected time series'
);
my $tempTS = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Temporary time series'
);
my $tempEventList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Temporary event list'
);
# my $lccorr_succeeded = doCommand(
# 'lccorr_pcms'
# , srctsset => $srcRateSet
# , eventset => $forTimeFilteredEvList
# , withbkgtsset => 'yes'
# , bkgtsset => $bgdRateSet
# , subtractbkg => 'yes'
# , treatvignet => 'no'
# , treatfiltertrans => 'no'
# , treatquantumeff => 'no'
# , treatcosmicrays => 'yes'
# , treatdeadtime => 'yes'
# , treatgti => 'yes'
# , treatalias => 'no'
# , applycorrections => 'yes'
# , outset => $bkgSubtractedRateSet
# , srcweightstyle => 'deduce' # this choice not supported in 'classical' lccorr.
# , srcra => $src_ra
# , srcdec => $src_dec
# , srcdeducestyle => 'events' # this parameter not supported in 'classical' lccorr.
#
# The following are parameters added to lccorr_pcms to compensate for it not reading the DSS of the rate file $srcRateSet. It is up to YOU to make sure that the parameter values correctly map the rate file event selections.
# , patternlo => 0
# , patternhi => $patternHi
# , flagmask => $eventFlagMask
# , pilos => 200
# , pihis => 12000
# , withgtis => 'yes'
# , gtitabsnode0 => [@gtiTabList]
# , rawylo => $rawylo
# , srcregiontab => "$srcRegionSet:REGION"
# , bkgregiontab => "$bkgRegionSet:REGION"
# , tempset => $tempTS
# , tempeventset => $tempEventList
#
# );
#
# replace lccorr_pcms by epiclccor
#############################################
# # x2, x4 REBINNING --------------------------------------
# # iterate over x2, x4 binning factors
#
# foreach my $rebin (qw(2 4)) {
#
# info("Generating x${rebin} rebinned source timeseries...");
#
# my $srcRateSetN = newFile(
# class => 'intermediate'
# , instrument => thisInstrument
# , src_num => $id
# , exp_id => $exp_id
# , content => "Uncorrected source time series (x${rebin} rebin)"
# );
#
# my $bgdRateSetN = newFile(
# class => 'intermediate'
# , instrument => thisInstrument
# , src_num => $id
# , exp_id => $exp_id
# , content => "Uncorrected background time series (x${rebin} rebin)"
# );
#
# doCommand(
# 'evselect'
# , table => "${forTimeFilteredEvList}:EVENTS"
# , writedss => 'yes'
# , updateexposure => 'yes'
# , expression => "region($srcRegionSet,X,Y)"
# , energycolumn => 'PI'
# , withrateset => 'yes'
# , rateset => $srcRateSetN
# , timebinsize => $rebin * $$timeBinSizeCol[$i]
# , maketimecolumn => 'yes'
# , makeratecolumn => 'yes'
# )
# or return exception();
#
# doCommand(
# 'evselect'
# , table => "${forTimeFilteredEvList}:EVENTS"
# , writedss => 'yes'
# , updateexposure => 'yes'
# , expression => "region($bkgRegionSet,X,Y)"
# , energycolumn => 'PI'
# , withrateset => 'yes'
# , rateset => $bgdRateSetN
# , timebinsize => $rebin * $$timeBinSizeCol[$i]
# , maketimecolumn => 'yes'
# , makeratecolumn => 'yes'
# )
# or return exception();
#
# # check whether output lightcurve has sufficient bins, skip if not
#
# my $num_bins = numberFITSRows(file => $srcRateSetN
# , extension => 'RATE');
#
# unless($num_bins >= 5)
# {
# info("x${rebin} rebinned lightcurve contains insufficient bins ($num_bins), skipping...");
# next;
# }
#
# info("Generating x${rebin} rebinned lightcurves...");
#
# my $bkgSubtractedRateSetN = newFile(
# class => 'intermediate'
# , instrument => thisInstrument
# , src_num => $id
# , exp_id => $exp_id
# , content => "Exposure-corrected time series (x${rebin} rebin)"
# );
#
# my $tempTSN = newFile(
# class => 'intermediate'
# , instrument => thisInstrument
# , src_num => $id
# , exp_id => $exp_id
# , content => "Temporary time series (x${rebin} rebin)"
# );
#
# my $tempEventListN = newFile(
# class => 'intermediate'
# , instrument => thisInstrument
# , src_num => $id
# , exp_id => $exp_id
# , content => "Temporary event list (x${rebin} rebin)"
# );
#
## my $lccorr_succeededN = doCommand(
## 'lccorr_pcms'
## , srctsset => $srcRateSetN
## , eventset => $forTimeFilteredEvList
## , withbkgtsset => 'yes'
## , bkgtsset => $bgdRateSetN
## , subtractbkg => 'yes'
## , treatvignet => 'no'
## , treatfiltertrans => 'no'
## , treatquantumeff => 'no'
## , treatcosmicrays => 'yes'
## , treatdeadtime => 'yes'
## , treatgti => 'yes'
## , treatalias => 'no'
## , applycorrections => 'yes'
## , outset => $bkgSubtractedRateSetN
## , srcweightstyle => 'deduce' # this choice not supported in 'classical' lccorr.
## , srcra => $src_ra
## , srcdec => $src_dec
## , srcdeducestyle => 'events' # this parameter not supported in 'classical' lccorr.
##
### The following are parameters added to lccorr_pcms to compensate for it not reading the DSS of the rate file $srcRateSet. It is up to YOU to make sure that the parameter values correctly map the rate file event selections.
## , patternlo => 0
## , patternhi => $patternHi
## , flagmask => $eventFlagMask
## , pilos => 200
## , pihis => 12000
## , withgtis => 'yes'
## , gtitabsnode0 => [@gtiTabList]
## , rawylo => $rawylo
## , srcregiontab => "$srcRegionSet:REGION"
## , bkgregiontab => "$bkgRegionSet:REGION"
## , tempset => $tempTSN
## , tempeventset => $tempEventListN
## );
##
## replace lccorr_pcms by epiclccor
#
# my $lccorr_succeededN = doCommand(
# 'epiclccorr'
# , srctslist => $srcRateSetN
# , eventlist => $forTimeFilteredEvList
# , withbkgset => 'yes'
# , bkgtslist => $bgdRateSetN
# , applyabsolutecorrections => 'no'
# , outset => $bkgSubtractedRateSetN
# );
#
# };
# # END REBINNING --------------------------------------
##############################################
# LC Corrections
# check results of x1 binning, proceed
my $process_ts_further = 0; # default # if $process_ts_further = 0 => epiclccorr does NOT work => NOT process further
if ( !readFITSKeyword( file => $srcRateSet, extension => 'RATE', keyword => 'EXPOSURE' )
|| !readFITSKeyword( file => $bgdRateSet, extension => 'RATE', keyword => 'EXPOSURE' ) ) {
info("No EXPOSURE keyword found in extension RATE of $srcRateSet or $bgdRateSet. Skipping further processing for this source.");
}
else
{
my $lccorr_succeeded = doCommand(
'epiclccorr'
, srctslist => $srcRateSet
, eventlist => $forTimeFilteredEvList
, withbkgset => 'yes'
, bkgtslist => $bgdRateSet
, applyabsolutecorrections => 'no'
, outset => $bkgSubtractedRateSet
);
if (!$lccorr_succeeded || !fileExists(file => $bkgSubtractedRateSet))
{
info("Lccorr failed for source $id");
# Skip rest of TS processing for this source, go straight to
# start making spectral products.
}
else
{
# lccorr succeeded, process further:
# I don't think this is needed anymore.
# doCommand(
# 'evselect'
# , table => $bkgSubtractedRateSet
# , keepfilteroutput => 'yes'
# , withfilteredset => 'yes'
# , filteredset => $tempTS
# , writedss => 'no'
# , updateexposure => 'no'
#FISH - spong evselect and n columns test below
# , expression => '!(isNull(RATE) && RATE > 0)'
# ) or return exception();
# Make sure the time series has non-null/zero bins before continuing.
my @rate = readFITSColumn(file => $bkgSubtractedRateSet
, extension => 'RATE'
, column => 'RATE'
);
if (grep { $_ } @rate)
{
$process_ts_further = 1;
}
else
{
info("Time series has no rows with non-null values of RATE. Not processed further.");
}
}
}
if ($process_ts_further) # if $process_ts_further = 1 => epiclccorr worked => process further
{
# gd - to add src flags:
doCommand( 'addattribute'
,attributename => [qw(SRCNUM SRC_RA SRC_DEC SRC_CRA SRC_CDEC EP_FLAG PN_FLAG M1_FLAG M2_FLAG EXTRRAD)]
,attributetype => [qw(integer real real real real string string string string real)]
,attributecomment => [
'EPIC source list number'
,'RA of source'
,'DEC of source'
,'Corrected RA of source'
,'Corrected DEC of source'
,'EP automatic source flag'
,'PN automatic source flag'
,'M1 automatic source flag'
,'M2 automatic source flag'
,'Source extraction radius (0.05 arcsec)'
]
,integervalue => $id
,realvalue => [ $$srcRaCol[$i] , $$srcDecCol[$i] , $$srcCorrRaCol[$i] , $$srcCorrDecCol[$i] , $source_radius[$id] ]
,stringvalue => [ $$epFlagCol[$i] , $$pnFlagCol[$i], $$m1FlagCol[$i], $$m2FlagCol[$i] ]
, set => "$bkgSubtractedRateSet:RATE"
) or exception();
my $mergedGtis = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Merged GTIs'
);
my @gtiTableList = ("${bkgSubtractedRateSet}:SRC_GTIS"
, "${bkgSubtractedRateSet}:BKG_GTIS");
if ($flareGtiExists) {push @gtiTableList, "${flareGtifile}:STDGTI";}
doCommand(
'gtimerge'
, tables => [@gtiTableList]
, withgtitable => 'yes'
, gtitable => $mergedGtis
, mergemode => 'and'
)
or return exception();
$alignedGtis = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Aligned merged GTIs'
);
doCommand(
'gtialign'
, style => 'generic'
, ingtitable => "$mergedGtis:STDGTI"
, outgtitable => "$alignedGtis:STDGTI"
, tstable => "${bkgSubtractedRateSet}:RATE"
)
or return exception();
my $tempCorrectedRateSet = newFile(
class => 'intermediate'
, instrument => thisInstrument
, band => 8
, src_num => $id
, exp_id => $exp_id
, content => 'Source timeseries'
);
# DLG manual add of variability headers for Mantis #18 ----
# Cannot set null values, throws AttNameValueMismatch
# Set -1.0 default as clearly unphysical
doCommand( 'addattribute'
, attributename => [qw(AVRATE CHISQUAR CHI2PROB N_POINTS)]
, attributetype => [qw(real real real integer)]
, attributecomment => [
'Average rate in light curve'
, 'Chi squared value'
, 'Probability that CHISQUAR arises by chance'
, 'The number of data points'
]
, realvalue => [ -1.0, -1.0, -1.0 ]
, integervalue => [ 0 ]
, set => "$bkgSubtractedRateSet:RATE"
) or exception();
# ---------------------------------------------------------
# Call ekstest twice (see Mantis 85)
# chisqr should run on unbackgroud-subtracted light curve
# fvar should run on background-subtracted light curve
my $ekstest_output = doCommand(
'ekstest'
, set => $bkgSubtractedRateSet
, gtis => 'yes'
, gtiset => $alignedGtis
, chisquaretest => 'yes'
, fracvartest => 'no'
, screen => 'yes'
, newoutset => 'yes'
, outset => $tempCorrectedRateSet
)
or return exception();
my $ekstest_output = doCommand(
'ekstest'
, set => $tempCorrectedRateSet
, gtis => 'yes'
, gtiset => $alignedGtis
, chisquaretest => 'no'
, fracvartest => 'yes'
, screen => 'yes'
, newoutset => 'yes'
, outset => $tempCorrectedRateSet
)
or return exception();
# my $fftPsFile = newFile(
# class => 'intermediate'
# , instrument => thisInstrument
# , exp_id => $exp_id
# , src_num => $id
# , content => 'Time series FFT plot'
# , format => 'PS'
# );
#
#
# my $efftplot_succeeded = doCommand(
# 'efftplot'
# , infile => $tempCorrectedRateSet
# , xscale => 'log'
# , yscale => 'lin'
# , gtis => 'yes'
# , gtiset => $alignedGtis
# , plotdevice => '/vcps'
# , outfile => $fftPsFile
# )
# or return exception();
#
# if (!$efftplot_succeeded || !fileExists(file => $fftPsFile))
# {
# info("efftplot failed for source $id");
# }
} # end 'if lccorr worked' sequelae
# Now create spectral products:
my $srcSelExpression = "(region($srcRegionSet,X,Y))";
my $bkgSelExpression = "(region($bkgRegionSet,X,Y))";
if ($flareGtiExists) {
$srcSelExpression .= " && (GTI(${flareGtifile}:STDGTI,TIME))";
$bkgSelExpression .= " && (GTI(${flareGtifile}:STDGTI,TIME))";
}
my $tempSrcSpectrum = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Source spectrum'
);
my $tempBkgSpectrum = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id # newFile takes care of hexification
, exp_id => $exp_id
, content => 'Source background spectrum'
);
my $srcSpectrum = newFile(
class => 'product'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'EPIC source spectrum'
) or exception();
my $bkgSpectrum = newFile(
class => 'product'
, instrument => thisInstrument
, src_num => $id # newFile takes care of hexification
, exp_id => $exp_id
, content => 'EPIC source background spectrum'
);
my $arf = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id # newFile takes care of hexification
, exp_id => $exp_id
, content => 'EPIC ancillary response function'
);
my $arfProd = newFile(
class => 'product'
, instrument => thisInstrument
, src_num => $id # newFile takes care of hexification
, exp_id => $exp_id
, content => 'EPIC ancillary response function'
);
doCommand(
'especget'
, table => $forSpecFilteredEvList
, srcexp => $srcSelExpression
, backexp => $bkgSelExpression
, withfilestem => 'no'
, srcspecset => $tempSrcSpectrum
, bckspecset => $tempBkgSpectrum
, srcarfset => $arf
, witharfset => 'yes'
, withrmfset => 'yes' # but leave --rmfset at default (yeesh,
# what an interface!) - this makes the
# task look for a canned rmf.
, useodfatt => 'no'
, withbadpixcorr => 'yes'
, extendedsource => 'no'
)
or return exception();
# If something went wrong and $tempSrcSpectrum was not created, skip
# to the next source.
if ( !fileExists( file => $tempSrcSpectrum ) ) {
info("$tempSrcSpectrum not created... skipping to next source.");
next;
}
# Enter source-specific keywords into spectral files:
# add autoSrcFlags here as above - gd
doCommand( 'addattribute'
,attributename => [qw(SRCNUM SRC_RA SRC_DEC SRC_CRA SRC_CDEC EP_FLAG PN_FLAG M1_FLAG M2_FLAG EXTRRAD)]
,attributetype => [qw(integer real real real real string string string string real)]
,attributecomment => [
'EPIC source list number'
,'RA of source'
,'DEC of source'
,'Corrected RA of source'
,'Corrected DEC of source'
,'EP automatic source flag'
,'PN automatic source flag'
,'M1 automatic source flag'
,'M2 automatic source flag'
,'Source extraction radius (0.05 arcsec)'
]
,integervalue => $id
,realvalue => [ $$srcRaCol[$i] , $$srcDecCol[$i] , $$srcCorrRaCol[$i] , $$srcCorrDecCol[$i] , $source_radius[$id] ]
,stringvalue => [ $$epFlagCol[$i] , $$pnFlagCol[$i] , $$m1FlagCol[$i], $$m2FlagCol[$i] ]
, set => ["$tempSrcSpectrum:SPECTRUM" , "$tempBkgSpectrum:SPECTRUM", "$arf:SPECRESP" ]
) or exception();
# For unexplicable reasons, especget is capable of producing a spectrum with no
# EXPOSURE keyword under some circumstances (namely after the event which
# damaged MOS1 on 11 DEC 2012). This check is a bandaid -- if no EXPOSURE
# keyword is found, skip this source, since otherwise it can cause
# a pipeline exception after grppha which seems to create an uninitialized
# EXPOSURE value. Note that we copy the temporary spectrum to the final
# product spectrum in this case since normally it is the output of grppha
# which we are now skipping (the main thing grppha does is add some extra
# binning information).
if ( !readFITSKeyword( file => $tempSrcSpectrum, extension => 2, keyword => "EXPOSURE" ) ) {
info("No EXPOSURE keyword found in $tempSrcSpectrum. Skipping further processing for this source.");
# PPR#7345: Intermediate spectra with EXPOSURE absent or set to zero should not be copied to final products.
#
# "and copying to product spectrum $srcSpectrum");
# copyFile(source => "$tempSrcSpectrum", destination => $srcSpectrum);
next;
}
# especget in this configuration doesn't either create or
# need an RMF file: all it does is decide which of the canned
# RMF matrices at the soc (sort of a pseudo-ccf) is the
# correct one to use for this source, filemode etc. Especget
# writes the name of this file (just the file name, no path
# or url information) to the keyword RESPFILE in the source
# spectrum dataset header. Apparently this keyword is where
# xspec expects to find the RMF file name.
#
# This whole section is a bit ludicrous. All that is required
# is to make a ps plot of the spectrum. There was a sas task,
# especplot, designed to do this. However there seemed to be
# eternal problems with especplot. For political reasons it
# was deemed inadvisable for an alternate developer to code
# this task. So it was decided to use xspec to make this plot.
# However xspec requires an RMF. Hence it is necessary to keep
# a cache of ALL the canned rmfs - just to make a plot!!
#
# In order for the whole Heath Robinson contraption to work,
# the simple file name in the RESPFILE keyword must be
# overwritten by a string of the form private_pcms_rmf_path/
# filename. The the simple file name must be put back into
# RESPFILE, for the benefit of a later interactive user of
# the spectrum. Eesh what a horrid system.
# For each of the three files ARF, RMF and bkgSpectrum, we need to know both the path name (to pass to xspec, so it can find the files) and the bare file name (to write to the output spectrum). The ARF and bkgSpectrum path names we know already, they are contained respectively in the variables $arf and $bkgSpectrum. These have to be cut up to yield the bare file name. For the RMF on the other hand we can read the bare file name from the SPECTRUM header of $tempSrcSpectrum, where especget writes it. There is a special pcms routine findRmfFile to transform this into a full path name.
my $cannedRmfFileName = readFITSKeyword(
file => $tempSrcSpectrum
, extension => 'SPECTRUM'
, keyword => 'RESPFILE'
);
info("Canned RMF is $cannedRmfFileName");
# Find the local path name $cannedRmfPathName:
my $cannedRmfPathName = findRmfFile( file => $cannedRmfFileName );
unless ( $cannedRmfPathName )
{
exception('MISSINGRMF','Unable to locate requested rmf : ', $cannedRmfFileName );
}
info("Absolute path to canned RMF is $cannedRmfPathName");
# Before running xspec we need to make sure make sure that
# $cannedRmfPathName exists and is readable. Otherwise xspec
# will end in an endless loop. - DJHF This is done by findRmfFile()
# ie. it returns undef when it cannot find a FILE.
# Cut off the directory names from the ARF and bkgSpectrum files:
my @chunks;
@chunks = split('\/', $arf);
my $arfFileName = $chunks[$#chunks];
@chunks = split('\/', $bkgSpectrum);
my $bkgSpectrumFileName = $chunks[$#chunks];
## my $grpphaCommExpr;
## if (thisInstrument eq 'epn')
## {
## $grpphaCommExpr = "bad 0-63&bad 1575-1645&group min "
## ."20&chkey RESPFILE $cannedRmfPathName"."&exit";
## } elsif (thisInstrument eq 'emos1' || thisInstrument eq 'emos2')
## {
## $grpphaCommExpr = "bad 0-23&group min "
## ."20&chkey RESPFILE $cannedRmfPathName"."&exit";
## }
##########################
# Define $spcgrpCommExpr (setbad) for 'specgroup'
my $spcgrpCommExpr;
if (thisInstrument eq 'epn')
{
$spcgrpCommExpr = '0.0:0.32'; # Spec range around 8 KeV ( 7.875:8.230 ) removed as bad range
} elsif (thisInstrument eq 'emos1' || thisInstrument eq 'emos2')
{
$spcgrpCommExpr = '0.0:0.36';
}
##########################
# Hope I've at last got this right! IMS 20051114
# Try rebinning in place as temp file - DLG
## doCommand(
## 'grppha'
## , infile => $tempSrcSpectrum
## , outfile => "\!$tempSrcSpectrum"
## , chatter => 0
## , comm => $grpphaCommExpr
## )
## or return exception();
##########################
# Replace of 'grppha' by 'specgroup'. JVP 23-06-2014
#
doCommand(
'specgroup'
, spectrumset => $tempSrcSpectrum # ok. -> intermediate/ Source spectrum
, groupedset => "\!$tempSrcSpectrum" # ok. Overwrite the input
, overwrite => 'Y' # ok
, backgndset => $tempBkgSpectrum # ok
, withbgdset => 'Y' # ok
, mincounts => 20 # ok. 25 for pnTiming.
, withCounts => 'Y' # ok
, oversample => 3 # ok
, withoversampling => 'Y' # ok
, rmfset => $cannedRmfPathName # ok
, withrmfset => 'Y' # ok
, arfset => $arf # ok
, witharfset => 'Y' # ok
, setbad => $spcgrpCommExpr # ok Defined by energy ranges.
, units => 'KEV' # ok
)
or return exception();
##########################
# Calculate total good counts (QUALITY=0) in spectrum
# Doing this the hard way because FTOOLS won't do what I want ;p
my $specQuality = readFITSColumn( file => $tempSrcSpectrum
, extension => 'SPECTRUM'
, column => 'QUALITY'
);
my $specCounts = readFITSColumn( file => $tempSrcSpectrum
, extension => 'SPECTRUM'
, column => 'COUNTS'
);
my $specNumRows = numberFITSRows( file => $tempSrcSpectrum
, extension => 'SPECTRUM'
);
my $specTotalCounts = 0;
for ($i = 0; $i < $specNumRows; $i++) {
if ( $$specQuality[$i] == 0 ) {
$specTotalCounts += $$specCounts[$i];
}
}
info("Source $id: Got $specTotalCounts total good counts at QUALITY==0");
# create spectrum product filenames
# generate xspec filename xspecIIIIEEEESSS.ps
my $xspecname = sprintf( '%s%s%s%03X' ,'xspec' , thisInstrument , $exp_id , $id);
my $spectrumPsFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'Spectrum plot'
, format => 'PS'
, name => "$xspecname.ps" # XSPEC doesn't like long filenames
);
my $spectrumPdfFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id # newFile takes care of hexification
, content => 'EPIC source spectrum plot'
, format => 'PDF'
);
my $xpecCommandFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'XSPEC command file'
, format => 'ASCII'
);
if ( $specTotalCounts < 100 ) {
info("Source $id: Fewer than 100 good counts in source spectrum, skipping product generation");
} else {
info("Source $id: More than 100 good counts in source spectrum, generating products");
# copy intermediate ARF to product
copyFile(
source => $arf
, destination => $arfProd
);
# copy intermediate spectrum to product
copyFile(
source => $tempSrcSpectrum
, destination => $srcSpectrum
);
# copy intermediate bkg spectrum to product
copyFile(
source => $tempBkgSpectrum
, destination => $bkgSpectrum
);
# copy intermediate timeseries to product
my $tempCorrectedRateSet = findFile(
class => 'intermediate'
, instrument => thisInstrument
, band => 8
, src_num => $id
, exp_id => $exp_id
, content => 'Source timeseries'
);
my $correctedRateSet = newFile(
class => 'product'
, instrument => thisInstrument
, band => 8
, src_num => $id
, exp_id => $exp_id
, content => 'EPIC source timeseries'
);
if (fileExists(file => $tempCorrectedRateSet)) {
copyFile(
source => $tempCorrectedRateSet
, destination => $correctedRateSet
);
# copy source and background region files to product
if (fileExists( file => $srcRegionSet ) ) {
doCommand(
'ftappend'
, infile => "${srcRegionSet}[REGION]"
, outfile => $correctedRateSet
) or return exception();
doCommand(
'fparkey'
, value => "REGION_SRC"
, fitsfile => "${correctedRateSet}[REGION]"
, keyword => "EXTNAME"
) or return exception();
};
if (fileExists( file => $bkgRegionSet) ) {
doCommand(
'ftappend'
, infile => "${bkgRegionSet}[REGION]"
, outfile => $correctedRateSet
) or return exception();
doCommand(
'fparkey'
, value => "REGION_BKG"
, fitsfile => "${correctedRateSet}[REGION]"
, keyword => "EXTNAME"
) or return exception();
};
# generate timeseries PDF product
my $tsPsFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'Source timeseries plot'
, format => 'PS'
);
doCommand(
'elcplot'
, set => $correctedRateSet
, usegtiset => 'yes'
, useflareset => 'no'
, gtiset => $alignedGtis
, offset => 'yes'
, offsetstyle => 'auto'
, forcestart => 'yes'
, binsize => 1
, plotdevice => '/vcps'
, plotfile => $tsPsFile
)
or return exception();
my $tsPdfFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, band => 8
, src_num => $id # newFile takes care of hexification
, content => 'EPIC source timeseries plot'
, format => 'PDF'
);
PStoPDF(
source => $tsPsFile
, destination => $tsPdfFile
) or return exception();
# Make the DS9 region product file:
my $srcASCRegionFile = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'src ds9 region'
, format => 'ASCII'
);
my $bkgASCRegionFile = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'bkg ds9 region'
, format => 'ASCII'
);
my $ds9RegionFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'EPIC source ds9 region'
, format => 'ASCII'
);
catASCIIFiles(
infiles => [$srcASCRegionFile, $bkgASCRegionFile, $allSrcRegionFile]
, outfile => $ds9RegionFile
) or return exception();
#################
#
# generate FFT Timeseries products
#
my $fftPsFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'Time series FFT plot'
, format => 'PS'
);
my $efftplot_succeeded = doCommand(
'efftplot'
, infile => $correctedRateSet
, xscale => 'log'
, yscale => 'lin'
, gtis => 'yes'
, gtiset => $alignedGtis
, plotdevice => '/vcps'
, outfile => $fftPsFile
)
or return exception();
if (!$efftplot_succeeded || !fileExists(file => $fftPsFile))
{
info("efftplot failed for source $id");
}
# generate FFT PDF plot product
#my $fftPsFile = findFile(
# class => 'intermediate'
# , instrument => thisInstrument
# , exp_id => $exp_id
# , src_num => $id
# , content => 'Time series FFT plot'
# , format => 'PS'
# );
my $fftPdfFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id # newFile takes care of hexification
, content => 'EPIC source FFT plot'
, format => 'PDF'
);
if(fileExists(file => $fftPsFile)) {
PStoPDF(
source => $fftPsFile
, destination => $fftPdfFile
) or return exception();
}
#################
}; # end of fileExists( file => $tempCorrectedRateSet )
# Data source name without leading directories for use in spectrum title.
my $fname = (split /\//,$srcSpectrum)[-1];
$fname =~ s/\.FIT$|\.FTZ$//g;
my $title = "Source $id minus background ($fname)\n";
open(XSPLIS, ">$xpecCommandFile");
print XSPLIS "data $srcSpectrum\n";
print XSPLIS "ignore **-0.2 12.0-**\n";
print XSPLIS "ignore bad\n";
print XSPLIS "setplot en\n";
print XSPLIS "cpd /NULL\n";
print XSPLIS "setplot command log y on\n";
print XSPLIS "setplot command r y\n";
# print XSPLIS "setplot command ld\n";
print XSPLIS "setplot command r x 0.2 12\n";
print XSPLIS "setplot command lwidth 6\n";
print XSPLIS "setplot command la FILE \n";
print XSPLIS "setplot command la T $title \n";
print XSPLIS "setplot command hardcopy $spectrumPsFile/ps\n";
print XSPLIS "plot\n";
print XSPLIS "exit\n";
print XSPLIS "Y\n";
close(XSPLIS);
doCommand( "xspec - '$xpecCommandFile'")
or return exception();
#################################
#################################
# If the spectrum file SPECTRUM extension has an ANCRFILe keyword
# rewrite the FIT to FTZ to maintain post-compression usability.
#my $ancrfile = readFITSKeyword(
# file => $srcSpectrum
# , extension => 'SPECTRUM'
# , keyword => 'ANCRFILE'
# );
#if ( $ancrfile )
#{
# # Make the keyword refer to a compressed file as this is what the user will
# # ultimately have.
#
# my $ancrFileName = $arfProd;
# $ancrFileName =~ s/\.FIT/\.FTZ/g;
#
# doCommand('addattribute'
# ,set => "$srcSpectrum:SPECTRUM"
# ,attributename => 'ANCRFILE'
# ,attributetype => 'string'
# ,stringvalue => $ancrFileName
# ,attributecomment => "\"\\\"ancillary response\\\"\""
# ) or exception();
#}
# Now rewrite the $cannedRmfFileName back into the RESPFILE keyword!!!!
# doCommand(
# 'addattribute'
# , set => "${srcSpectrum}:PRIMARY"
##### whoops - addattribute is not yet able to write kwds to an extension header! 'PRIMARY' must be wrong in any case. But it seems not to matter here...
# , attributename => 'RESPFILE'
# , attributetype => 'string'
# , stringvalue => $cannedRmfFileName
# )
# or return exception();
# INSERTED BY DLG 26.01.09 ---------------------
# Manual modification of ANCRFILE and FITSFILE headers
# to reflect compressed FITS filenames
# $arfFileName =~ s/\.FIT/\.FTZ/g;
#my $ancrFileName = $arfProd;
#$ancrFileName =~ s/\.FIT/\.FTZ/g;
#$ancrFileName =~ s/product\///g;
#$bkgSpectrumFileName =~ s/\.FIT/\.FTZ/g;
#
#$grpphaCommExpr = "chkey ANCRFILE $ancrFileName"."&chkey RESPFILE "
# ."$cannedRmfFileName"."&chkey BACKFILE $bkgSpectrumFileName"."&exit";
# INSERTED BY DLG 26.01.09 ---------------------
#doCommand(
# 'grppha'
# , infile => $srcSpectrum
# , outfile => "\!$srcSpectrum"
# , chatter => 0
# , comm => $grpphaCommExpr
# )
# or return exception();
##########################
# Replace grppha by addattribute on changing filenames keywords
#
my $ancrFileKey = readFITSKeyword(
file => $srcSpectrum
, extension => 'SPECTRUM'
, keyword => 'ANCRFILE'
);
my $rmfFileKey = readFITSKeyword(
file => $srcSpectrum
, extension => 'SPECTRUM'
, keyword => 'RESPFILE'
);
my $backFileKey = readFITSKeyword(
file => $srcSpectrum
, extension => 'SPECTRUM'
, keyword => 'BACKFILE'
);
# Make the keyword refer to a compressed file as this is what the user will ultimately have.
my $ancrFileName = $arfProd;
$ancrFileName =~ s/\.FIT/\.FTZ/g;
$ancrFileName =~ s/product\///g;
my $rmfFileName = $cannedRmfFileName;
$bkgSpectrumFileName =~ s/\.FIT/\.FTZ/g;
if ( $ancrFileKey && $rmfFileKey && $backFileKey ) {
doCommand('addattribute'
,set => "$srcSpectrum:SPECTRUM"
,attributename => [qw(ANCRFILE RESPFILE BACKFILE)]
,attributetype => [qw(string string string)]
,stringvalue => [ $ancrFileName, $rmfFileName, $bkgSpectrumFileName ]
,attributecomment => [
'ancillary response'
, 'response matrix'
, 'background spectrum'
]
) or exception();
}
##########################
#################################
#################################
if (fileExists(file => $spectrumPsFile))
{
PStoPDF(
source => $spectrumPsFile
, destination => $spectrumPdfFile
)
or return exception();
}
else
{
info("Can't find $spectrumPsFile - so can't make PDF.");
}
} # end generating products for $specTotalCounts > 100
# if (!(fileExists(file => $spectrumPsFile) || $process_ts_further))
# {
# info("There is neither a PS spectrum product nor a time series with any");
# info("useful data; thus DS9-region product for this source won't be made.");
# next;
# }
} # end loop over sources
info("********** End of forked child process $$");
}; # End of the eval
# This is the end of the forked child process. Write the
# seqdb to a temporary file which we will use to
# communicate information back to the parent process
# (e.g., the names of new files that were created). If
# an exception was detected inside the eval {} block,
# we also record it to the seqdb file.
my $exception;
if ($@) {
$exception = $@;
}
writeChildSeqdb(exception => $exception);
# Return status of child process reflects occurence of an exception
if ($exception) {
exit -1;
} else {
exit 0;
}
} else {
# Couldn't fork
return exception("Couldn't fork: $!\n");
}
} # end loop over jobs
# Now wait for all the child jobs to finish.
# Check to see if any gave a bad exit status before returning.
$badexitstatus = waitChild( \@children );
# Necessary so that DB connections may be destroyed again by the parent
Pcmsdb::set_InactiveDestroy(0);
# Merge sequence information from the children
my $extraarrayref = mergeChildSeqdb(\@children);
if ($badexitstatus) {
# mergeChildSeqdb should have copied any exception in a child
# process into the parent Seqdb which will be picked up in the
# following exception() call.
return exception();
}
#-------------------------------------------------------------------------------------------------------------------
#
# Flags SPECTRA and TSERIES columns in final SourceList only if at least one epic camera
# has produced Spectrum + TS
#
# Find final Products files (Spectra and TS) for each EPIC instrument: emos1, emos2, epn
my @finalSpectra_files = findFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'EPIC source spectrum'
, 'format' => 'FITS'
);
my (@epn_spectra, @emos1_spectra, @emos2_spectra, @epic_spectra);
my $inst = thisInstrument;
# If there are no Spectra_files
if ($#finalSpectra_files<0) { info("DEBUG SRC_NUM No source spectra for $inst , exp $exp_id"); }
#
my %src_radius; # Hash to keep the pair SRCNUM => Exatraction region radius
foreach my $f (@finalSpectra_files){
next unless fileExists (file => $f);
my $src_num = readFITSKeyword(
file => $f
, extension => 'SPECTRUM'
, keyword => 'SRCNUM'
) or return exception();
#info("DEBUG Source: $src_num SPECTRA for $inst : $f");
if ($inst eq 'epn'){ push (@epn_spectra , $src_num);}
if ($inst eq 'emos1'){ push (@emos1_spectra , $src_num);}
if ($inst eq 'emos2'){ push (@emos2_spectra , $src_num);}
###################################
# Reading Source Extractoin region raius from Spectrum REGION extension
my $src_radius = readFITSColumn(
file => $f
, extension => "REGION"
, column => "R"
);
info("DEBUG : Source Extraction radius from Spectrum $f : $$src_radius[0]");
my @arr = ($inst, $exp_id, $$src_radius[0]);
$src_radius{$src_num} = \@arr;
### $src_radius{$src_num} = $$src_radius[0];
###################################
#
# Create EPIC thumbnail images with Src + Bkg extraction regions (implotregions)
#
my $bkgregfile = findFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $src_num
, content => 'EPIC source background spectrum'
, format => 'FITS'
);
#my $outThumbnail = newFile(
# class => 'intermediate'
# , instrument => thisInstrument
# , exp_id => $exp_id
# , src_num => $src_num
# , content => 'EPIC thumbnail image'
# , format => 'PNG'
# );
my $outThumbnail = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $src_num
, content => 'EPIC source spectrum'
, format => 'PNG'
);
if ( fileExists( file => $f ) && fileExists( file => $bkgregfile ) ){
doCommand('implotregions'
,imageset => $epicImage
,srcregfile => $f
,bkgregfile => $bkgregfile
,outputfile => $outThumbnail
#) or exception();
) or info("DEBUG : Some problem in 'implotregions'. EPIC thumbnail image not created for SRCNUM: $src_num");
} else {
info("DEBUG : Src spec $f or Bkg spec $bkgregfile not found. EPIC thumbnail image not created.");
}
###################################
}
#
# Intermediate ASCII file to store src_num ids of sources with Spectrum+TS for each epic Instrument and Exposure
my $epic_spectra_reg;
$epic_spectra_reg = findFile(
class => 'intermediate'
, instrument => 'epic'
, content => 'SRC_NUM with Spectra and TS per instrument'
, format => 'ASCII'
);
unless ( $epic_spectra_reg && fileExists( file => $epic_spectra_reg ) ){
info("DEBUG Creating file $epic_spectra_reg . Inst $inst , exp $exp_id");
$epic_spectra_reg = newFile(
class => 'intermediate'
, instrument => 'epic'
, content => 'SRC_NUM with Spectra and TS per instrument'
, format => 'ASCII'
);
}
# Write info in intermediate file:
my @lines = ();
if ($inst eq 'epn'){
info("DEBUG SRC_NUM ids for sources with epn spectrum: @epn_spectra");
push (@lines, "SRC_NUM ids for sources with epn spectrum Exp. $exp_id : @epn_spectra\n"); }
if ($inst eq 'emos1'){
info("DEBUG SRC_NUM ids for sources with emos1 spectrum: @emos1_spectra");
push (@lines, "SRC_NUM ids for sources with emos1 spectrum Exp. $exp_id : @emos1_spectra\n"); }
if ($inst eq 'emos2'){
info("DEBUG SRC_NUM ids for sources with emos2 spectrum: @emos2_spectra");
push (@lines, "SRC_NUM ids for sources with emos2 spectrum Exp. $exp_id : @emos2_spectra\n"); }
#push (@lines, "SRC_NUM ids for sources with EPIC spectrum Exp. $exp_id : @epic_spectra\n");
writeASCIIFile( name => $epic_spectra_reg
, text => \@lines
, append => 1
);
# Read and get final SRC_NUM list from intermediate ASCII file:
my @records = readASCIIFile(name => $epic_spectra_reg);
my @arr;
foreach my $line (@records){
my @fields = split / /, (split /: /, $line)[1];
push (@arr, @fields);
}
# @epic_spectra uniq sorted:
@epic_spectra = sort { $a <=> $b } ( do { my %seen; grep { !$seen{$_}++ } @arr } );
info("DEBUG SRC_NUM ids for sources with at least one spectrum in any EPIC camera: @epic_spectra");
# Write results to hash %src_spectra_tseries
foreach my $source ( sort { $a <=> $b } keys(%src_spectra_tseries) ){
if ( $source ~~ @epic_spectra ){ # If $id source is in @selectedSources list => True, True
##@{$src_spectra_tseries{$source}} = (1,1); # ( SPECTRA, TSERIES ) = ( True, True )
#info("Source SRC_NUM = $source flagged to SPECTRA, TSERIES = True, True , because src has final products");
#if ( $src_radius{$source} ){
## @{$src_spectra_tseries{$source}} = (1,1, $src_radius{$source} );
##@{$src_spectra_tseries{$source}} = (1,1, @{$src_radius{$source}}[0], @{$src_radius{$source}}[1], @{$src_radius{$source}}[2] );
# Check if there is TimeSeries available from any EPIC camera
my $src_tseries_emos1 = findFile(
class => 'product'
, instrument => 'emos1'
, src_num => $source
, content => 'EPIC source timeseries'
, 'format' => 'FITS'
);
my $src_tseries_emos2 = findFile(
class => 'product'
, instrument => 'emos2'
, src_num => $source
, content => 'EPIC source timeseries'
, 'format' => 'FITS'
);
my $src_tseries_epn = findFile(
class => 'product'
, instrument => 'epn'
, src_num => $source
, content => 'EPIC source timeseries'
, 'format' => 'FITS'
);
if ( fileExists(file => $src_tseries_emos1) || fileExists(file => $src_tseries_emos2) || fileExists(file => $src_tseries_epn)){
##info("DEBUG ==> Time Series exists: $src_tseries_emos1, $src_tseries_emos2, $src_tseries_epn");
@{$src_spectra_tseries{$source}} = (1,1, $inst, $exp_id, @{$src_radius{$source}}[2] );
} else {
@{$src_spectra_tseries{$source}} = (1,0, $inst, $exp_id, @{$src_radius{$source}}[2] );
}
#}
} else { # other sources NOT in @selectedSources => False, False
@{$src_spectra_tseries{$source}} = (0,0, $inst, $exp_id, 0); # ( SPECTRA, TSERIES ) = ( False, False )
#info("Source SRC_NUM = $source flagged to SPECTRA, TSERIES = False, False , because src has not final products");
}
}
## DEBUG =================================
foreach my $k ( sort { $a <=> $b } keys(%src_spectra_tseries) ){
info("DEBUG ==> SRCNUM: $k RADIUS: @{$src_spectra_tseries{$k}}");
}
## DEBUG =================================
# Write final flags : SPECTRA and TSERIES from hash %src_spectra_tseries to columns in source list: $summaryList
my $summaryList = findFile(
class => "product"
, instrument => 'epic'
, content => 'EPIC summary source list'
, format => 'FITS'
) or return exception();
# Modify Source List
my $summaryList_modification = writeFITSfromHash($summaryList , %src_spectra_tseries);
if ( $summaryList_modification ){
info("writeFITS: SPECTRA and TSERIES flags written to Source List.");
}
####################################################
# Create intermediate summaryList to write the Src extraction radius for implot:
my $summaryListIntermediate = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'EPIC summary source list for implot'
#, format => 'ASCII'
);
copyFile(
source => $summaryList
, destination => $summaryListIntermediate
) or return exception();
# Modify Intermediate Source List for implot
my $summaryList_modification = writeRadiusfromHash($summaryListIntermediate , %src_spectra_tseries);
####################################################
#-------------------------------------------------------------------------------------------------------------------
return success();
}
# -----------------------------------------------------------------------------
# Use SAS task ecoordconv to convert RA/Dec to detector coordinates
# Task returns WALLOFTEXT rather than sensible parameters, so this will be ugly hack :p
#
sub coordConv
{
# coordConv: convert RA/Dec to detector coordinates
# coordConv(image, ra, decl)
#
#
info("--- Coordinate conversion subroutine ---");
my ($image, $ra, $decl) = @_;
# check input validity
if ($image eq '') { Exception->error('Module', "coordConv: Image not defined") };
if ($ra !~ /^\d+(\.\d*)?/) { Exception->error('Module', "coordConv: RA $ra not valid") };
if ($decl !~ /^-?\d+(\.\d*)?/) { Exception->error('Module', "coordConv: Dec $decl not valid") };
# invoke ecoordconv and capture output
my $cmd = "ecoordconv imageset=$image x=$ra y=$decl coordtype=eqpos withcoords=yes";
info("CMD: $cmd");
my $output = `$cmd`;
info("$output");
# Search output for XY values (yuk)
my ($x,$y);
($x, $y) = ($output =~ /X:\s+Y:\s+([\d.]+)\s+([\d.]+)/);
my $sccd;
( $sccd ) = ($output =~ /centred on CCD: (\d+)/);
# info("coordConv returns source CCD: $sccd");
# Return reference to X, Y values
return ($x, $y, $sccd);
}
# -----------------------------------------------------------------------------
# Add row to region file. Original region file is created from scratch
#
sub writeRegionRow2
{
my ($file, $shape, $x, $y, $r_i, $r_o) = @_;
my ($data, $ff, $nr, $status);
my $exp_id = exposureID( instrument => thisInstrument
, stream => thisStream
);
# Generate data template for REGION columns
# Populate with data
my @shape = ( $shape );
my @x = ( $x );
my @y = ( $y );
my @r;
if ( !defined $r_o || $r_o == 0 ){
@r = ( $r_i );
} else {
@r = ( $r_i, $r_o );
}
my $dim_r = scalar @r;
info("DEBUG - Creating FITS REGION file for $shape - r: @r - Dimen r: $dim_r");
# Check if region file exists
# if file doesn't exist, create it
if (!fileExists( file => $file )) {
info("writeRegionRow(): Region file does not exist: creating it");
# Create a new REGION FITS file
fits_create_file($ff, $file, $status);
# Empty PRIMARY extension
fits_create_img($ff, 16, 0, 0, $status);
# Includes CREATOR and DATE in the fits file
fits_write_key($ff, TSTRING, 'CREATOR', 'EPICSourceProducts.pm', 'name of code which created this file', $status);
fits_write_date($ff, $status); # current system date as a character string in 'yyyy-mm-ddThh:mm:ss' format
# Create empty FITS table with the proper columns
if ( $shape =~ /ANNULUS/i) { # Shape = annulus -> 2 dimensions for R and 5 columns
fits_create_tbl($ff, BINARY_TBL, 0, 5, ['SHAPE','X','Y', 'R', 'COMPONENT'], ['16A', 'E', 'E', '2E', 'B'], ['', '', '', '', ''], 'REGION', $status);
FITSException->error($status,'Error creating ', file => $file) if ($status);
} else { # Shape != annulus -> 1 dimensions for R and 4 columns
fits_create_tbl($ff, BINARY_TBL, 0, 4, ['SHAPE','X','Y', 'R'], ['16A', 'E', 'E', 'E'], ['', '', '', ''], 'REGION', $status);
FITSException->error($status,'Error creating ', file => $file) if ($status);
}
$ff->close_file($status);
# write expected headers to REGION extension
# see ims_ascregion_mod.f90 in SAS 'region' source for details - DLG
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'HDUVERS'
, attributetype => 'string'
, stringvalue => '1.0.0'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'HDUCLASS'
, attributetype => 'string'
, stringvalue => 'ASC'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'HDUCLAS1'
, attributetype => 'string'
, stringvalue => 'REGION'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'HDUCLAS2'
, attributetype => 'string'
, stringvalue => 'STANDARD'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'MTYPE1'
, attributetype => 'string'
, stringvalue => 'pos'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'MFORM1'
, attributetype => 'string'
, stringvalue => 'X,Y'
) or exception();
}
# Open file for READWRITE
fits_open_file($ff, $file, READWRITE, $status);
FITSException->error($status,'Error reading ', file => $file) if ($status);
# Move to REGION extension
$ff->movnam_hdu(ANY_HDU,'REGION', 0, $status);
FITSException->error($status,'Error accessing extension ', file => $file, extension => 'REGION') if ($status);
# Get number of rows
$ff->get_num_rows($nr, $status);
FITSException->error($status,'Error reading num_rows ', file => $file, extension => 'REGION') if ($status);
info("DEBUG: REGION file $file: writing $shape, $x, $y, $r_i, $r_o firstrow $nr");
# write data to table
# advance one row from $nr
$ff->write_col(TSTRING, 1, $nr+1, 1, 1, \@shape, $status);
FITSException->error($status,'Error writing SHAPE ', file => $file, extension => 'REGION') if ($status);
$ff->write_col(TFLOAT, 2, $nr+1, 1, 1, \@x, $status);
FITSException->error($status,'Error writing X ', file => $file, extension => 'REGION') if ($status);
$ff->write_col(TFLOAT, 3, $nr+1, 1, 1, \@y, $status);
FITSException->error($status,'Error writing Y ', file => $file, extension => 'REGION') if ($status);
$ff->write_col(TFLOAT, 4, $nr+1, 1, $dim_r, \@r, $status); # R can be 2 dimensions array if shape = annulus
FITSException->error($status,'Error writing R array ', file => $file, extension => 'REGION') if ($status);
# # If shape = !CIRCLE and $file has 5 columns
# if ( $shape =~ /!CIRCLE/ && hasFITSKeyword( file => $file, extension => 'REGION', keyword => 'TTYPE5') ) {
# my @compo = ( 1 );
# If $file has 5 columns ==> All elements (ANNULUS region and !CIRCLE's) belong to the same COMPONENT = 1
if ( hasFITSKeyword( file => $file, extension => 'REGION', keyword => 'TTYPE5') ) {
my @compo = ( 1 );
$ff->write_col(TBYTE, 5, $nr+1, 1, 1, \@compo, $status);
FITSException->error($status,'Error writing Component column', file => $file, extension => 'REGION') if ($status);
}
# Close file
$ff->close_file($status);
}
# -----------------------------------------------------------------------------
# Add row to region file. Original region file from template 'extraction region file' which is emptied before populated
#
sub writeRegionRow
{
my ($file, $shape, $x, $y, $r) = @_;
my ($data, $ff, $nr, $status);
my $exp_id = exposureID( instrument => thisInstrument
, stream => thisStream
);
# Generate data template for REGION columns
# Populate with data
my @shape = ( $shape );
my @x = ( $x );
my @y = ( $y );
my @r = ( $r );
# Check if region file exists
# if file doesn't exist, create it
if (!fileExists( file => $file )) {
info("writeRegionRow(): Region file does not exist: creating it");
# find master region file
my $extrRegion = findFile(
class => "intermediate"
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'extraction region file'
, format => 'FITS'
) or return exception();
# dummy fselect to generate a file with empty REGION table
doCommand(
'fselect'
, infile => $extrRegion
, outfile => $file
, expr => "SHAPE .eq. 'FOO'"
);
# write expected headers to REGION extension
# see ims_ascregion_mod.f90 in SAS 'region' source for details - DLG
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'HDUVERS'
, attributetype => 'string'
, stringvalue => '1.0.0'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'HDUCLASS'
, attributetype => 'string'
, stringvalue => 'ASC'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'HDUCLAS1'
, attributetype => 'string'
, stringvalue => 'REGION'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'HDUCLAS2'
, attributetype => 'string'
, stringvalue => 'STANDARD'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'MTYPE1'
, attributetype => 'string'
, stringvalue => 'pos'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'MFORM1'
, attributetype => 'string'
, stringvalue => 'X,Y'
) or exception();
}
# Open file for READWRITE
fits_open_file($ff, $file, READWRITE, $status);
FITSException->error($status,'Error reading ', file => $file) if ($status);
# Move to REGION extension
$ff->movnam_hdu(ANY_HDU,'REGION', 0, $status);
FITSException->error($status,'Error accessing extension ', file => $file, extension => 'REGION') if ($status);
# Get number of rows
$ff->get_num_rows($nr, $status);
FITSException->error($status,'Error reading num_rows ', file => $file, extension => 'REGION') if ($status);
info("DEBUG: REGION file $file: writing $shape, $x, $y, $r firstrow $nr");
# write data to table
# advance one row from $nr
$ff->write_col(TSTRING, 1, $nr+1, 1, 1, \@shape, $status);
FITSException->error($status,'Error writing SHAPE ', file => $file, extension => 'REGION') if ($status);
$ff->write_col(TFLOAT, 2, $nr+1, 1, 1, \@x, $status);
FITSException->error($status,'Error writing X ', file => $file, extension => 'REGION') if ($status);
$ff->write_col(TFLOAT, 3, $nr+1, 1, 1, \@y, $status);
FITSException->error($status,'Error writing Y ', file => $file, extension => 'REGION') if ($status);
$ff->write_col(TFLOAT, 4, $nr+1, 1, 1, \@r, $status);
FITSException->error($status,'Error writing R ', file => $file, extension => 'REGION') if ($status);
# Close file
$ff->close_file($status);
}
# -----------------------------------------------------------------------------
sub writeASCIIRegion
{
# write DS9 ASCII region file
#
# circular or annulus region in detector coordinates
# writeAsciiRegion($file, $x, $y, $r, $text, $append)
my ( $name, $shape, $x, $y, $r_i, $r_o, $colour, $label, $append ) = @_;
if ($name eq '') { Exception -> error( 'Module', "writeAsciiRegion: filename not defined") };
if ($x !~ /^\d+(\.\d*)?/) { Exception -> error( 'Module', "writeAsciiRegion: detector X $x not valid") };
if ($y !~ /^\d+(\.\d*)?/) { Exception -> error( 'Module', "writeAsciiRegion: detector Y $y not valid") };
if ($r_i !~ /^\d+(\.\d*)?/) { Exception -> error( 'Module', "writeAsciiRegion: region radius R_i $r_i not valid") };
if ($r_o !~ /^\d+(\.\d*)?/) { Exception -> error( 'Module', "writeAsciiRegion: region radius R_o $r_o not valid") };
# Generate parameter strings
my $labstring = "";
if ($label ne "") { $labstring = " text={${label}}"; }
if ($colour eq '') { $colour = 'white'; }
my $colstring = "# color=".$colour;
# Generate array for text output
my @t = ();
if ($append == 0) {
push (@t, "# Region file format: DS9 version 3.0\n");
push (@t, "global font = \"helvetica 10 normal\" edit=0 move=0 delete=1 include=1 fixed=0\n");
if ( $shape =~ /CIRCLE/i) {
push (@t, "detector;circle($x, $y, $r_i) ${colstring}${labstring}\n");
} elsif ( $shape =~ /ANNULUS/i) {
push (@t, "detector;annulus($x, $y, $r_i, $r_o) ${colstring}${labstring}\n");
}
### push (@t, "detector;circle($x, $y, $r) ${colstring}${labstring}\n");
} else {
push (@t, "detector;-circle($x, $y, $r_i) ${colstring}${labstring}\n");
}
# Dump results to ASCII file
writeASCIIFile( name => $name
, text => \@t
, append => $append
);
my $rtn = 1;
return $rtn;
}
# -----------------------------------------------------------------------------
sub readFITStoHash {
my @f = @_;
my (@src, @spectra, @tseries);
my %out_hash;
my $nsrc=numberFITSRows(file => @f,
extension => "SRCLIST");
my $SRC_NUM_ = readFITSColumn(
file => @f
, extension => "SRCLIST"
, column => "SRC_NUM"
);
my $SPECTRA_ = readFITSColumn(
file => @f
, extension => "SRCLIST"
, column => "SPECTRA"
);
my $TSERIES_ = readFITSColumn(
file => @f
, extension => "SRCLIST"
, column => "TSERIES"
);
my $j;
for ($j = 0; $j < $nsrc; $j++ ) {
my $src = $$SRC_NUM_[$j];
my $spectra = $$SPECTRA_[$j];
my $tseries = $$TSERIES_[$j];
my @arr = ($spectra, $tseries);
$out_hash{$src} = \@arr;
}
return %out_hash;
}
# -----------------------------------------------------------------------------
sub writeFITSfromHash {
my ($file, %d) = @_;
my (@src, @spectra, @tseries);
if (!fileExists( file => $file )) {
info("EPIC summary source list does not exist.");
}
if (!hasFITSExtension(file => $file, extension => "SRCLIST")){
info("File $file doesn't contain extension SRCLIST - Source list won't be modified.");
}
####################################################################
#
# Read srcList and get SRC_NUM, SPECTRA and TSERIES columns
my $nsrc=numberFITSRows(file => $file,
extension => "SRCLIST");
my $SRC_NUM_ = readFITSColumn(
file => $file
, extension => "SRCLIST"
, column => "SRC_NUM"
);
my $SPECTRA_ = readFITSColumn(
file => $file
, extension => "SRCLIST"
, column => "SPECTRA"
);
my $TSERIES_ = readFITSColumn(
file => $file
, extension => "SRCLIST"
, column => "TSERIES"
);
my $j;
for ($j = 0; $j < $nsrc; $j++ ) {
my $src = $$SRC_NUM_[$j];
my $spectra = $$SPECTRA_[$j];
my $tseries = $$TSERIES_[$j];
####################################################################
#
# Change the SPECTRA and TSERIES values to False if %d{$src} = (0,0)
#
if ( !@{$d{$src}}[0] && !@{$d{$src}}[1] ){ # spec=0, tseries=0, F,F
$spectra = $tseries = 0;
} elsif ( @{$d{$src}}[0] && !@{$d{$src}}[1] ){ # spec=1, tseries=0, T,F
$spectra = 1;
$tseries = 0;
} else {
$spectra = $tseries = 1; # spec=1, tseries=1, T,T
}
info("DEBUG SRC_NUM writeFITS: Source= $src . Changing SPECTRA= $spectra , TSERIES= $tseries"); ############################
####################################################################
push (@src , $src);
push (@spectra , $spectra);
push (@tseries , $tseries);
}
####################################################################
#
# Modify the srcList with the new SPECTRA and TSERIES column values
my @new_SPECTRA = @spectra;
my @new_TSERIES = @tseries;
my $data;
$$data{SPECTRA} = [@new_SPECTRA];
$$data{TSERIES} = [@new_TSERIES];
$$data{-column}{SPECTRA} = { 'width' => 4,
#'colnum' => 6,
'ttype' => 'SPECTRA',
'repeat' => 1,
'typecode' => 42
};
$$data{-column}{TSERIES} = { 'width' => 4,
#'colnum' => 7,
'ttype' => 'TSERIES',
'repeat' => 1,
'typecode' => 42
};
my $colname = $$data{colname};
push (@$colname, 'SPECTRA', 'TSERIES');
$$data{colname} = [ @$colname ];
writeFITSColumn(
file => $file
, extension => "SRCLIST"
, column => $data
);
}
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
sub writeRadiusfromHash {
my ($file, %d) = @_;
my (@src, @radius);
if (!fileExists( file => $file )) {
info("EPIC summary source list does not exist.");
}
if (!hasFITSExtension(file => $file, extension => "SRCLIST")){
info("File $file doesn't contain extension SRCLIST - Source list won't be modified.");
}
####################################################################
#
# Read srcList and get SRC_NUM
my $nsrc=numberFITSRows(file => $file,
extension => "SRCLIST");
my $SRC_NUM_ = readFITSColumn(
file => $file
, extension => "SRCLIST"
, column => "SRC_NUM"
);
my $j;
for ($j = 0; $j < $nsrc; $j++ ) {
my $src = $$SRC_NUM_[$j];
####################################################################
#
# Insert Column Radius to the Intermediate Source List
#
my $radium = @{$d{$src}}[4];
info("DEBUG SRC_NUM writeFITS: Source= $src . RADIUS= $radium"); ############################
####################################################################
push (@src , $src);
push (@radius , $radium);
}
## DEBUG ----------------------------------------------------
#info("DEBUG ==> Interm. Source list : $file");
#foreach my $k (@radius){
# info("DEBUG ==> RADIUS: $k");
#}
## DEBUG ----------------------------------------------------
####################################################################
#
# Modify the srcList with the new RADIUS column values
my @new_RADIUS = @radius;
my $data;
$$data{RADIUS} = [@new_RADIUS];
$$data{-column}{RADIUS} = { 'width' => 4,
#'colnum' => 7,
'ttype' => 'RADIUS',
'repeat' => 1,
'typecode' => 42
};
my $colname = $$data{colname};
push (@$colname, 'RADIUS');
$$data{colname} = [ @$colname ];
writeFITSColumn(
file => $file
, extension => "SRCLIST"
, column => $data
);
}
# -----------------------------------------------------------------------------
1;