Download

# 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.14';
$version = $VERSION;
$author  = 'Jose Vicente Perea';
$date    = '2016-12-12';

#
# ChangeLog
# =========
#
# 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 $totalcounts = readFITSKeyword(
                                   file => $forSpecFilteredEvList
                                   , extension => 'EVENTS'
                                   , keyword => "NAXIS2"
                                   );
				   
    my $timebinsize = $livetime/$totalcounts*100. ;

    if ($timebinsize < ($livetime/1000.)) {
	$timebinsize = $livetime/1000. ;
    }


    # 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])'
              , 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])'
              , withtimeranges => 'Y'
              , timemin => $mintime
              , timemax => $maxtime
              , timebinsize => $timebinsize
              , maketimecolumn => "yes"
              , makeratecolumn => "yes"
              )
        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;