# This module generates source products for pn timing-mode
# observations
package pnTimingProducts;
use strict;
use English;
use Carp;
use vars
qw(@ISA $VERSION $name $author $date $version @instList $numberOfStreams);
@ISA = qw(Module);
use ModuleResources;
# Declare identity, version, author, date, etc.
$name = __PACKAGE__;
$VERSION = '0.16';
$version = $VERSION;
$author = 'Jose Vicente Perea';
$date = '2018-05-21';
#
# ChangeLog
# =========
#
# Version 0.16 - 2018-05-21 (JVP)
# ------------
#
# + Time bin size for PN Light Curves fix to 20 times the Frametime: timebinsize = 20 * frmtime
#
# Version 0.15 - 2018-04-24 (JVP)
# ------------
#
# + Create Light Curves per Energy Band. Bands: 2, 3, 4 and 5 (Timing case)
# 'Uncorrected source/background time series per EnergyBand'
# Those Light Curves are then merged by 'srcTimeSeriesEnergyBands.py' from LCMerge.pm
#
# Version 0.14 - 2016-12-12 (JVP)
# ------------
#
# + Do not ignore ~ 8 KeV spectral region (7.0:9.0) for spectrum products
#
# Version 0.13 - 2016-02-29 (JVP)
# ------------
#
# + Special region filter function 'box' replaced by open RAWX regions in src and bkg expressions.
#
# Version 0.12 - 2014-11-24 (JVP)
# ------------
#
# + Calls to ( 'evselect' + 'backscale' + 'arfgen' ) replaced by only one call to 'especget'.
# As a consecuence the PN-CannedRMF is not harcoded, but it is found by 'especget'.
# especget (withrmfset=yes and rmfset not set) search for the PN-CannedRMF which best suits with the spectrum epoch.
#
# Version 0.11 - 2014-08-08 (JVP)
# ------------
#
# + Correction in subroutine 'coordConvRADec'. Declination number is +<>, -<> or <>.
#
# Version 0.10 - 2014-08-07 (JVP)
# ------------
#
# + Range of valid maximum peak changed. Valid range for RAWX '15 -- 52', instead '12 -- 52'.
#
#
# Version 0.09 - 2014-08-05 (JVP)
# ------------
#
# + Subroutine 'coordConv' changed to 'coordConvRADec' to avoid confusion with 'coordConv' in EPICSourceProducts.pm module.
#
# Version 0.08 - 2014-07-30 (JVP)
# ------------
#
# + Object coordinates in LC plots are obtained from the PN Timing collapsed image where
# RAWX = maximum in the histogram ($wmean), RAWY = 190 ; which is the source position.
#
# Version 0.07 - 2014-06-23 (JVP)
# ------------
#
# + Re-naming ANCRFILE, RESPFILE and BACKFILE keywords with only one 'addattribute' call
#
# Version 0.06 - 2014-06-17 (JVP,PR)
# ------------
#
# + Skip product generation on <100 good counts
# + RMF file is not calculated anymore. A canned version is taken instead:
# 'epn_ti40_sdY9.rmf' for TIME.EL mode and 'epn_bu23_sdY9.rmf' for BURST.EL
# + Set bad data in 0.-0.5 and 7.0-9.0 keV energy ranges
# + Re-naming ANCRFILE, RESPFILE and BACKFILE keywords in products:
# Absolute path removed and .FIT to .FTZ in file extension.
# + Object coordinates in LC plots are RA_NOM, DEC_NOM, but labeled as 'Source RA/DEC_corr'. To be corrected.
# + Background flaring filtering is not attempted in pn Timing modes. The reason is that a significant fraction of flares found
# in the background region (RAWX in [2,6]) at energies > 10keV are actually due to flares in the source brightness.
#
# Version 0.05 - 2014-05-30 (JVP,PR)
# ------------
#
# + New products creation: EPIC ancillary response function, EPIC source spectrum plot,
# EPIC source timeseries and timeseries plot and EPIC source FFT timeseries plot.
#
# Version 0.04 - 2014-05-27 (JVP,PR)
# ------------
#
# + Light Curves: 'Background subtracted Light Curve' to product. 'EPIC source timeseries and timeseries plot'. Mantis issue 7185.
#
# Version 0.03 - 2014-05-12 (JVP,PR)
# ------------
#
# + Histogram creation: Call to 'fhisto' replaced by 'evselect'.
# FITS file extension: HISTO. Column names: RAWX, COUNTS
# Light Curves creation: 'Background subtracted Light Curve' filename to intermediate.
# Mantis issue 7185.
#
# Version 0.02 - 2014-04-24 (JVP)
# ------------
#
# + Change EPIC source and background spectra filenames to product. Mantis issue 7185.
#
# Version 0.01 - 2014-03-11 (EC)
# ------------
#
# + Initial version based on script extract_timing.csh from MG
#
# Declare list of instruments this module is interested in
@instList = qw(epn);
# Number of streams
sub numberOfStreams
{
return numberOfExposures();
}
# Rules method
sub evaluateRules
{
#
return ignore()
if ignored(
module => 'pnEvents'
, instrument => thisInstrument
, stream => thisStream
);
# ignore for slew processing
return ignore() if ($ENV{'PCMS_ISSLEW'});
#
start()
if complete(
module => 'pnEvents'
, instrument => thisInstrument
, stream => thisStream
);
}
# Action method
sub performAction
{
info("Module version number: $version");
# Work out which exposure this stream maps on to
my $exp_id = exposureID(
instrument => thisInstrument
, stream => thisStream
);
info("Exposure ID: $exp_id");
#
my $filter = getExposureProperty(
exp_id => $exp_id
, name => 'filter'
);
my $mode = getExposureProperty(
exp_id => $exp_id
, name => 'mode'
);
info("Filter: $filter");
info("Instrument mode: $mode");
# Find event files for this exposure
my $evFile = findFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => '0'
, content => 'EPIC timing mode event list'
);
# Return immediately if there are no event lists
return success() if ! $evFile;
# nothing to do unless we're in timing mode
my $datamode = readFITSKeyword(
file => $evFile
, extension => 5
, keyword => "DATATYPE"
);
if ( ! ($datamode eq "TIME.EL" || $datamode eq "BURST.EL") ) {
info("DATATYPE in ${evFile}[5] is $datamode, not TIME.EL or BURST.EL. Skipping.");
return success();
}
# Determine the column with the greatest brightness
my $colHistFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
#, src_num => $srcnum
, content => "Column histogram"
);
# doCommand(
# 'fhisto'
# , infile => $evFile
# , outfile => $colHistFile
# , column => "RAWX"
# , binsz => "INDEF"
# ) or return exception();
doCommand(
'evselect'
, table => $evFile
, expression => "(FLAG==0)&&(PATTERN<=4)&&(PI IN [500:7000])&&(RAWX in [1:63])"
, withhistogramset => 'yes'
, histogramset => $colHistFile
, histogramcolumn => "RAWX"
, histogrambinsize => 1
, withhistoranges => "yes"
, histogrammin => 1
, histogrammax => 64
) or return exception();
my $hist = readFITSTable(
file => $colHistFile
, extension => 'HISTO'
, colname => [qw( RAWX COUNTS )]
);
if ( $#{$hist->{RAWX}} != 63 ) {
Exception->error('FHISTO',
"The last row index of $colHistFile should be 63 but it is $#{$hist->{X}}"
);
}
# Skip the first column, when searching for the peak. Also
# skip producing products if the peak is too near the edges.
my $wmean=0;
my $maxhist=0;
my $specTotalCounts = 0;
for my $i ( 1 .. 63 ) {
my $x = ${$hist->{RAWX}}[$i];
my $y = ${$hist->{COUNTS}}[$i];
if ($y > $maxhist) {
$maxhist = $y;
$wmean = $x;
}
# Total counts from histogram
$specTotalCounts += $y;
}
info("ExposureID $exp_id: Got $specTotalCounts total good counts at FLAG==0, PATTERN<=4, etc");
# Return immediately if total number of good counts < 100
my $countsThreshold = 100;
if ( $specTotalCounts < $countsThreshold ) {
info("ExposureID $exp_id: Got $specTotalCounts total good counts which is < $countsThreshold. Skipping.");
return success();
}
# Return immediately if Peak RAWX occurs outside of range 15 -- 52
if ( ($wmean > 52) || ($wmean < 15) ) {
info("Peak RAWX occurs at $wmean which is not in the range 15 -- 52. Skipping.");
return success();
}
#***************************************************************************************************
#***************************************************************************************************
## info("Extracting source spectrum centered on RAWX=$wmean");
##
## my $forSpecFilteredEvList = newFile(
## class => 'intermediate'
## , instrument => thisInstrument
## , exp_id => $exp_id
## #, src_num => $srcnum
## , content => 'Event list for spectra'
## );
##
## doCommand(
## 'evselect'
## , withfilteredset => 'Y'
## , keepfilteroutput => 'Y'
## , table => "$evFile"
## , filteredset => $forSpecFilteredEvList
## , expression => "((RAWX,RAWY) IN box($wmean,100,9,100,0))"
## )
## or return exception();
##
##
## # The (non-binned)source spectrum. Should be a product in the future?
## my $srcSpecFile = newFile(
## class => 'intermediate'
## , instrument => thisInstrument
## , exp_id => $exp_id
## #, src_num => $srcnum
## , content => 'PN TIMING SOURCE SPECTRUM'
## );
##
## doCommand(
## 'evselect'
## , table => $forSpecFilteredEvList
## , withspectrumset => 'Y'
## , spectrumset => $srcSpecFile
## , energycolumn => 'PI'
## , filtertype => 'expression'
## , expression => '(PATTERN<=4 && FLAG==0)'
## , withspecranges => 'Y'
## , specchannelmin => 0
## , specchannelmax => 20479
## , spectralbinsize => 5
## )
## or return exception();
##
## doCommand(
## 'backscale'
## , spectrumset => $srcSpecFile
## , badpixlocation => $evFile
## )
## or return exception();
##
## info("Extracting the background spectrum ...");
##
## my $forBgSpecFilteredEvList = newFile(
## class => 'intermediate'
## , instrument => thisInstrument
## , exp_id => $exp_id
## #, src_num => $srcnum
## , content => 'Event list for BG spectra'
## );
##
## doCommand(
## 'evselect'
## , withfilteredset => 'Y'
## , keepfilteroutput => 'Y'
## , table => "$evFile"
## , filteredset => $forBgSpecFilteredEvList
## , expression => "((RAWX,RAWY) IN box(4,100,2,100,0))"
## )
## or return exception();
##
##
## # The background spectrum.
## my $bgSpecFile = newFile(
## class => 'product'
## , instrument => thisInstrument
## , exp_id => $exp_id
## , src_num => '0'
## , content => 'EPIC SOURCE BACKGROUND SPECTRUM'
## );
##
## doCommand(
## 'evselect'
## , table => $forBgSpecFilteredEvList
## , withspectrumset => 'Y'
## , spectrumset => $bgSpecFile
## , energycolumn => 'PI'
## , filtertype => 'expression'
## , expression => '(PATTERN<=4 && FLAG==0)'
## , withspecranges => 'Y'
## , specchannelmin => 0
## , specchannelmax => 20479
## , spectralbinsize => 5
## )
## or return exception();
##
## doCommand(
## 'backscale'
## , spectrumset => $bgSpecFile
## , badpixlocation => $evFile
## )
## or return exception();
##
##
##
##
## ########################
## #info("Calculating RMF ...");
## #
## #my $rmfFile = newFile(
## # class => 'intermediate'
## # , instrument => thisInstrument
## # , exp_id => $exp_id
## # #, src_num => $srcnum
## # , content => 'RMF for spectrum'
## # );
## #
## #doCommand(
## # 'rmfgen'
## # , rmfset => $rmfFile
## # , spectrumset => $srcSpecFile
## # )
## # or return exception();
## #
## ########################
##
## # Canned RMF files for PN Timing mode. Hard-coded at 12-06-2014 (JVP)
## # Timing = "TIME.EL" or Burst = "BURST.EL"
## my $cannedRmfFileName;
## if ( $datamode eq "TIME.EL" ) {
## $cannedRmfFileName = 'epn_ti40_sdY9.rmf';
## } elsif ( $datamode eq "BURST.EL" ) {
## $cannedRmfFileName = 'epn_bu23_sdY9.rmf';
## }
##
## info("Obtaining path for canned RMF file $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");
##
## my $rmfFile = $cannedRmfPathName ;
##
## ########################
##
## info("Calculating ARF ...");
##
## # The response matrix.
## my $ResponseMatrix = newFile(
## class => 'product'
## , instrument => thisInstrument
## , exp_id => $exp_id
## , src_num => '0'
## , content => "EPIC ANCILLARY RESPONSE FUNCTION"
## );
##
## doCommand(
## 'arfgen'
## , spectrumset => $srcSpecFile
## , arfset => $ResponseMatrix
## , withrmfset => 'Y'
## , rmfset => $rmfFile
## , badpixlocation => $evFile
## , detmaptype => 'flat'
## )
## or return exception();
#***************************************************************************************************
#***************************************************************************************************
# 'especget' is replacing ( 2 x evselect + backscale + arfgen ) for SourceSpec + BkgSpec
#
info("Extracting source spectrum centered on RAWX=$wmean and the background spectrum ...");
# The (non-binned) source spectrum. Should be a product in the future?
my $srcSpecFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'PN TIMING SOURCE SPECTRUM'
);
# The background spectrum.
my $bgSpecFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => '0'
, content => 'EPIC SOURCE BACKGROUND SPECTRUM'
);
# The response matrix.
my $ResponseMatrix = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => '0'
, content => "EPIC ANCILLARY RESPONSE FUNCTION"
);
doCommand(
'especget'
, table=> $evFile
#, srcexp => "(((RAWX,RAWY) IN box($wmean,100,9,100,0)) && PATTERN<=4 && FLAG==0)"
#, backexp => "(((RAWX,RAWY) IN box(4,100,2,100,0)) && PATTERN<=4 && FLAG==0)"
, srcexp => "(( (RAWX >= ($wmean-9)) && (RAWX <= ($wmean+9)) ) && PATTERN<=4 && FLAG==0)" # 'box' replaced by open RAWX range
, backexp => "(( (RAWX >= (4-2)) && (RAWX <= (4+2)) ) && PATTERN<=4 && FLAG==0)" # 'box' replaced by open RAWX range
, withfilestem => 'N'
, srcspecset => $srcSpecFile
, bckspecset => $bgSpecFile
, srcarfset => $ResponseMatrix
, withrmfset => 'Y'
)
or return exception();
# especget (withrmfset=yes and rmfset not set) search for the CannedRMF which best suits with the spectrum epoch
# and it is written in the RESPFILE keyword of the spectrum file.
# After that the CannedRMF file name is read from that keyword from $srcSpecFile
my $cannedRmfFileName = readFITSKeyword(
file => $srcSpecFile
, 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");
my $rmfFile = $cannedRmfPathName ;
#***************************************************************************************************
#***************************************************************************************************
info("Spectral rebinning ...");
# The binned source spectrum.
my $binnedSrcSpecFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => '0'
, content => 'EPIC SOURCE SPECTRUM'
);
doCommand(
'specgroup'
, spectrumset => $srcSpecFile
, groupedset => $binnedSrcSpecFile
, backgndset => $bgSpecFile
, withbgdset => 'Y'
, mincounts => 25
, withCounts => 'Y'
, oversample => 3
, withoversampling => 'Y'
, rmfset => $rmfFile
, withrmfset => 'Y'
, arfset => $ResponseMatrix
, witharfset => 'Y'
#, setbad => '0:0.5,7.0:9.0'
, setbad => '0:0.5' # Spec range around 8 KeV ( 7.0:9.0 ) removed as bad range
, units => 'KEV'
)
or return exception();
# Spectrum Plot - start - ###########################################################################################
# generate xspec filename xspecIIIIEEEESSS.ps
#my $xspecname = sprintf( '%s%s%s%03X' ,'xspec' , thisInstrument , $exp_id , $id);
my $xspecname = sprintf( '%s%s%s%03X' ,'xspec' , thisInstrument , $exp_id , '0' );
my $spectrumPsFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
#, src_num => $srcnum
, 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 => '0' # 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 => $srcnum
, content => 'XSPEC command file'
, format => 'ASCII'
);
# Input spectrum : $binnedSrcSpecFile from 'specgroup' task
# Data source name without leading directories for use in spectrum title.
my $fname = (split /\//,$binnedSrcSpecFile)[-1];
$fname =~ s/\.FIT$|\.FTZ$//g;
my $title = "Source 0 minus background ($fname)\n"; # Source = 0 <= there is not source in timing mode
# Bad channels defined in previous specgroup call
open(XSPLIS, ">$xpecCommandFile");
print XSPLIS "data $binnedSrcSpecFile\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 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();
############################################################################
# Re-naming ANCRFILE, RESPFILE and BACKFILE in spectrum KEYWORDS
# - Path to files removed in order to work in the cwd
# - FIT to FTZ in filenames
### ANCRFILE
# Cut off the directory names from the ARF and bkgSpectrum files:
my @chunks;
@chunks = split('\/', $ResponseMatrix);
my $arfFileName = $chunks[$#chunks];
@chunks = split('\/', $bgSpecFile);
my $backFileName = $chunks[$#chunks];
my $ancrFileKey = readFITSKeyword(
file => $binnedSrcSpecFile
, extension => 'SPECTRUM'
, keyword => 'ANCRFILE'
);
### RESPFILE
my $rmfFileKey = readFITSKeyword(
file => $binnedSrcSpecFile
, extension => 'SPECTRUM'
, keyword => 'RESPFILE'
);
### BACKFILE
my $backFileKey = readFITSKeyword(
file => $binnedSrcSpecFile
, extension => 'SPECTRUM'
, keyword => 'BACKFILE'
);
#####################
if ( $ancrFileKey && $rmfFileKey && $backFileKey ) {
my $ancrFileName = $arfFileName;
$ancrFileName =~ s/\.FIT/\.FTZ/g;
my $rmfFileName = $cannedRmfFileName;
$backFileName =~ s/\.FIT/\.FTZ/g;
doCommand('addattribute'
,set => "$binnedSrcSpecFile:SPECTRUM"
,attributename => [qw(ANCRFILE RESPFILE BACKFILE)]
,attributetype => [qw(string string string)]
,stringvalue => [ $ancrFileName, $rmfFileName, $backFileName ]
,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.");
}
# Spectrum Plot - end - ###########################################################################################
# Light curves (timeseries) - start - ####################################################################################
my $forSpecFilteredEvList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Event list for spectra'
);
doCommand(
'evselect'
, withfilteredset => 'Y'
, keepfilteroutput => 'Y'
, table => "$evFile"
, filteredset => $forSpecFilteredEvList
#, expression => "((RAWX,RAWY) IN box($wmean,100,9,100,0))"
, expression => "( (RAWX >= ($wmean-9)) && (RAWX <= ($wmean+9)) )" # 'box' replaced by open RAWX range
)
or return exception();
my $forBgSpecFilteredEvList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Event list for BG spectra'
);
my $srcRateFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
#, src_num => $srcnum
, content => 'PN TIMING SOURCE LIGHT CURVE'
);
my $bgRateFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
#, src_num => $srcnum
, content => 'PN TIMING BACKGROUND LIGHT CURVE'
);
my $livetime = readFITSKeyword(
file => $forSpecFilteredEvList
, extension => 'EVENTS'
, keyword => "LIVETIME"
);
my $mintime = readFITSKeyword(
file => $forSpecFilteredEvList
, extension => 'EVENTS'
, keyword => "TSTART"
);
my $maxtime = readFITSKeyword(
file => $forSpecFilteredEvList
, extension => 'EVENTS'
, keyword => "TSTOP"
);
my $frmtime = readFITSKeyword(
file => $forSpecFilteredEvList
, extension => 'EVENTS'
, keyword => "FRMTIME"
);
my $totalcounts = readFITSKeyword(
file => $forSpecFilteredEvList
, extension => 'EVENTS'
, keyword => "NAXIS2"
);
# my $timebinsize = $livetime/$totalcounts*100. ;
#
# if ($timebinsize < ($livetime/1000.)) {
# $timebinsize = $livetime/1000. ;
# }
# ------------ Time Bin size for Light Curves:
# timebinsize = 20 * FRMTIME PN
#
$frmtime = $frmtime * 1e-03; # [seconds]. Original FITS keyword in miliseconds
my $timebinsize = 20 * $frmtime ;
info("Source 0 time bin size for LC. Instr.: epn timing, Frametime: $frmtime, TimeBinSize = $timebinsize");
# Could the following 2 x evselect calls be reduced to only 1 call in the future ?
doCommand(
'evselect'
, table => $forSpecFilteredEvList
, withrateset => 'Y'
, rateset => $srcRateFile
, timecolumn => 'TIME'
, filtertype => 'expression'
, expression => '(PATTERN<=4)&&(FLAG==0)&&(PI in [500:7000])'
, energycolumn => 'PI'
, withtimeranges => 'Y'
, timemin => $mintime
, timemax => $maxtime
, timebinsize => $timebinsize
, maketimecolumn => "yes"
, makeratecolumn => "yes"
)
or return exception();
# Could the following 2 x evselect calls be reduced to only 1 call in the future ?
doCommand(
'evselect'
, withfilteredset => 'Y'
, keepfilteroutput => 'Y'
, table => "$evFile"
, filteredset => $forBgSpecFilteredEvList
#, expression => "((RAWX,RAWY) IN box(4,100,2,100,0))"
, expression => "( (RAWX >= (4-2)) && (RAWX <= (4+2)) )" # 'box' replaced by open RAWX range
)
or return exception();
doCommand(
'evselect'
, table => $forBgSpecFilteredEvList
, withrateset => 'Y'
, rateset => $bgRateFile
, timecolumn => 'TIME'
, filtertype => 'expression'
, expression => '(PATTERN<=4)&&(FLAG==0)&&(PI in [500:7000])'
, energycolumn => 'PI'
, withtimeranges => 'Y'
, timemin => $mintime
, timemax => $maxtime
, timebinsize => $timebinsize
, maketimecolumn => "yes"
, makeratecolumn => "yes"
)
or return exception();
############################################################################
#
# Time Series per Energy Band :
#
# PN Timing all Energy range : [500:7000] =>
# Band 1 is excluded and upper limit of Band 5 is changed to 7000
#
info("Time Series per Energy Bands -- SRC = 0 (timing), EXPOSURE = $exp_id");
foreach my $band ( 2, 3, 4, 5 ){
my ( $pilo, $pihi );
my $srcRateSetEband = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => '0'
, exp_id => $exp_id
, band => $band
, content => 'Uncorrected source time series per EnergyBand'
);
my $bgdRateSetEband = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => '0'
, exp_id => $exp_id
, band => $band
, content => 'Uncorrected background time series per EnergyBand'
);
# EnergyBands ranges:
$pilo = energyBand( $band, 0 );
$pihi = energyBand( $band, 1 );
my $ID_BAND = $band;
if ( $band == 5 ){
$pihi = 7000 ; # Upper limit of Band 5 changed to 7000
$ID_BAND = '100' ; # ID_BAND changed to 100 because the upper limit has been modified
}
info("Band: $ID_BAND . pilo = $pilo, pihi = $pihi");
doCommand(
'evselect'
, table => $forSpecFilteredEvList
, withrateset => 'Y'
, rateset => $srcRateSetEband
, timecolumn => 'TIME'
, filtertype => 'expression'
, expression => "(PATTERN<=4)&&(FLAG==0) && (PI in ($pilo:$pihi])"
, energycolumn => 'PI'
, withtimeranges => 'Y'
, timemin => $mintime
, timemax => $maxtime
, timebinsize => $timebinsize
, maketimecolumn => "yes"
, makeratecolumn => "yes"
)
or return exception();
doCommand(
'evselect'
, table => $forBgSpecFilteredEvList
, withrateset => 'Y'
, rateset => $bgdRateSetEband
, timecolumn => 'TIME'
, filtertype => 'expression'
, expression => "(PATTERN<=4)&&(FLAG==0) && (PI in ($pilo:$pihi])"
, energycolumn => 'PI'
, withtimeranges => 'Y'
, timemin => $mintime
, timemax => $maxtime
, timebinsize => $timebinsize
, maketimecolumn => "yes"
, makeratecolumn => "yes"
)
or return exception();
# Add the ID_BAND to the headers
doCommand(
'addattribute', set => "$srcRateSetEband:RATE"
, attributename => 'ID_BAND'
, attributetype => 'integer'
, attributecomment => '"\"EPIC Energy Band range ID\""'
, integervalue => $ID_BAND
) or return exception();
doCommand(
'addattribute', set => "$bgdRateSetEband:RATE"
, attributename => 'ID_BAND'
, attributetype => 'integer'
, attributecomment => '"\"EPIC Energy Band range ID\""'
, integervalue => $ID_BAND
) or return exception();
}
############################################################################
# Light Curve background subtracted.
# EPIC source timeseries and EPIC source timeseries plot
my $bkgSubtractedRateFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => '0'
, content => 'EPIC source timeseries'
);
doCommand(
'epiclccorr'
, srctslist => $srcRateFile
, eventlist => $evFile
, withbkgset => 'yes'
, bkgtslist => $bgRateFile
, applyabsolutecorrections => 'no'
, outset => $bkgSubtractedRateFile
)
or return exception();
###################
# Add coordinates SRC_CRA,SRC_CDEC to $bkgSubtractedRateFile from PN timing collapsed image (evselect -> eccordconv)
#
my $ra_nom = readFITSKeyword(
file => $srcRateFile
, extension => 'PRIMARY'
, keyword => "RA_NOM"
);
my $dec_nom = readFITSKeyword(
file => $srcRateFile
, extension => 'PRIMARY'
, keyword => "DEC_NOM"
);
info("RA_NOM, DEC_NOM from $srcRateFile = $ra_nom , $dec_nom");
# Sky coordinates obtained from collapsed image ('evselect') by using 'ecoordconv'
#
my $imgFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'PN TIMING source image'
);
doCommand(
'evselect'
, table => $evFile
, withimageset => 'yes'
, imageset => $imgFile
, xcolumn => 'DETX'
, ycolumn => 'DETY'
, ximagebinsize => 20
, yimagebinsize => 20
, imagebinning => 'binSize'
)
or return exception();
#### 'ecoordconv' task returns WALLOFTEXT rather than sensible parameters (not confirmed yet),
#### so $src_cra,$src_cdec is captured from the standard output with sub 'coordConvRADec' (see below)
####
####my $src_cra;
####my $src_cdec;
####doCommand(
#### 'ecoordconv'
#### , imageset => $imgFile
#### , coordtype => 'raw'
#### , withcoords => 'true'
#### , x => $wmean
#### , y => 190
#### , withccd => 'true'
#### , ccdno => 4
#### , ra_out => $src_cra
#### , dec_out => $src_cdec
#### )
#### or return exception();
####
my ($src_cra, $src_cdec, $sccd) = coordConvRADec($imgFile, $wmean, 190);
info("SRC_CRA, SRC_CDEC from ecoordconv in $imgFile = $src_cra , $src_cdec");
if ( $src_cra && $src_cdec ) {
doCommand('addattribute'
,set => "$bkgSubtractedRateFile:RATE"
, attributename => [qw(SRC_CRA SRC_CDEC)]
, attributetype => [qw(real real)]
, realvalue => [ $src_cra, $src_cdec ]
, attributecomment => [
'[deg] RA in source position'
, '[deg] Dec in source position'
]
) or exception();
}
info("Coordinates RA,Dec from ecoordconv harcoded as SRC_CRA,SRC_CDEC in file = $bkgSubtractedRateFile");
###########################
#############################################
# Getting alignedGtis
#
my $mergedGtis = newFile(
class => 'intermediate'
, instrument => thisInstrument
#, src_num => $srcnum
, exp_id => $exp_id
, content => 'Merged GTIs'
);
my @gtiTableList = ("${bkgSubtractedRateFile}:SRC_GTIS"
, "${bkgSubtractedRateFile}:BKG_GTIS");
# Pending to un-comment if it is needed (JVP 2014-05-27)
## if ($flareGtiExists) {push @gtiTableList, "${flareGtifile}:STDGTI";}
doCommand(
'gtimerge'
, tables => [@gtiTableList]
, withgtitable => 'yes'
, gtitable => $mergedGtis
, mergemode => 'and'
)
or return exception();
my $alignedGtis = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
#, src_num => $srcnum
, content => 'Aligned merged GTIs'
);
doCommand(
'gtialign'
, style => 'generic'
, ingtitable => "$mergedGtis:STDGTI"
, outgtitable => "$alignedGtis:STDGTI"
, tstable => "${bkgSubtractedRateFile}:RATE"
)
or return exception();
#############################################
my $tsPsFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
#, src_num => $srcnum
, content => 'Source timeseries plot'
, format => 'PS'
);
doCommand(
'elcplot'
, set => $bkgSubtractedRateFile
, usegtiset => 'yes'
, useflareset => 'no'
, gtiset => $alignedGtis
, offset => 'yes'
, offsetstyle => 'auto'
, forcestart => 'yes'
, binsize => 1
, plotdevice => '/vcps'
, plotfile => $tsPsFile
)
or return exception();
# Light Curve plot. PDF format
# band = 0 in the first version. To be decided by PR
my $tsPdfFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, band => '0'
, src_num => '0'
, content => 'EPIC source timeseries plot'
, format => 'PDF'
);
PStoPDF(
source => $tsPsFile
, destination => $tsPdfFile
) or return exception();
##########################################
# Light Curve FFT Plot
my $tempCorrectedRateFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, band => 8
#, src_num => $srcnum
, 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 => "$bkgSubtractedRateFile: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 => $bkgSubtractedRateFile
, gtis => 'yes'
, gtiset => $alignedGtis
, chisquaretest => 'yes'
, fracvartest => 'no'
, screen => 'yes'
, newoutset => 'yes'
, outset => $tempCorrectedRateFile
)
or return exception();
my $ekstest_output = doCommand(
'ekstest'
, set => $tempCorrectedRateFile
, gtis => 'yes'
, gtiset => $alignedGtis
, chisquaretest => 'no'
, fracvartest => 'yes'
, screen => 'yes'
, newoutset => 'yes'
, outset => $bkgSubtractedRateFile # overwrite again $bkgSubtractedRateFile
)
or return exception();
# generate FFT timeseries PDF product
my $fftPsFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
#, src_num => $srcnum
, content => 'Time series FFT plot'
, format => 'PS'
);
my $efftplot_succeeded = doCommand(
'efftplot'
, infile => $bkgSubtractedRateFile
, xscale => 'log'
, yscale => 'lin'
, gtis => 'yes'
, gtiset => $alignedGtis
, plotdevice => '/vcps'
, outfile => $fftPsFile
);
if (!$efftplot_succeeded || !fileExists(file => $fftPsFile))
{
info("efftplot failed for source $exp_id");
}
my $fftPdfFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, band => '0'
, src_num => '0'
, content => 'EPIC SOURCE FFT PLOT'
, format => 'PDF'
);
PStoPDF(
source => $fftPsFile
, destination => $fftPdfFile
) or return exception();
# Light curves (timeseries) - end - ####################################################################################
# Everything OK
return success();
}
# -----------------------------------------------------------------------------
# Use SAS task ecoordconv to convert detector coordinates to RA/Dec
# Task returns WALLOFTEXT rather than sensible parameters, so this will be ugly hack :p
#
sub coordConvRADec
{
# coordConvRADec: convert detector coordinates X,Y to RA/Dec
# coordConvRADec(image, X, Y)
#
info("--- Coordinate conversion subroutine ---");
my ($image, $X, $Y) = @_;
# check input validity
if ($image eq '') { Exception->error('Module', "coordConvRADec: 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") };
if ( $X !~ /^[-]?\d+$/ ) { Exception->error('Module', "coordConvRADec: X $X not valid") };
if ( $Y !~ /^[-]?\d+$/ ) { Exception->error('Module', "coordConvRADec: Y $Y not valid") };
# invoke ecoordconv and capture output
my $cmd = "ecoordconv imageset=$image x=$X y=$Y coordtype=raw withcoords=yes withccd=yes ccdno=4";
info("CMD: $cmd");
my $output = `$cmd`;
info("$output");
# Search output for RA,Dec values
my ($RA,$Dec);
($RA, $Dec) = ($output =~ /RA:\s+DEC:\s+(\d+\.\d+)\s+((?:\+|-|)\d+\.\d+)/);
info("coordConvRADec returns RA: DEC: $RA $Dec");
my $sccd;
( $sccd ) = ($output =~ /centred on CCD: (\d+)/);
info("coordConvRADec returns source CCD: $sccd");
# Return reference to RA, Dec values
return ($RA, $Dec, $sccd);
}
# -----------------------------------------------------------------------------
1;