# 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;