package EPICSourceProducts; use strict; use English; use Carp; use Regexp::Common qw(boolean); use vars qw(@ISA $VERSION $name $author $date $version @instList $numberOfStreams); @ISA = qw(Module); use ModuleResources; use ModuleUtil; use Astro::FITS::CFITSIO qw(TRUE FALSE); # Declare identity, version, author, date, etc. $name = __PACKAGE__; $VERSION = '1.25'; $version = $VERSION; $author = 'Anja Schroeder,Ian Stewart,Duncan John Fyfe'; $date = '2006-10-13'; # # ChangeLog # ========= # # version 1.26 2009-01-26 DLG # ------------ # # + Corrected ANCRFILE and BACKFILE headers in output filename to reflect compressed FITS files # # # version 1.25 2006-10-13 DJF # ------------ # # + Make sure each image exists before we try to access fits keywords, # + Make sure we have some images to act on aswell. # # version 1.24 2006-09-20 DJF # ------------ # # + Record use of exposures for source specific product creation in the database. # # version 1.23 2006-09-15 DJF # ------------ # # + evaluateRules added test on ExpChooser. Exposures ignored by ExpChooser were being processed. This should not happen. # # version 1.22 - 2006-09-14 DJF # ------------ # # + Construct short unique filename to replace hardcoded xpsec.ps for xpsec postscript output. # + Apply isEmptyImage test to band 1-5 images. If any are empty do not proceed with source specific product creation. # # version 1.21 - 2006-09-05 DJF # ------------ # # + Use fmodhead to re-write ANCRFILE keyword in spectrum file SPECTRUM extension. the extension .FIT becomes .FTS to maintain post-compression usability. # # version 1.20 - 2006-08-08 ACS # ------------ # # + changed line width in xspec plot from 3 to 6 (for better conversion into gif-files) # # version 1.19 - 2006-07-13 DJF # ------------ # # + Removed # from source spectrum title - xpsec doesn't like it. # # version 1.18 - 2006-07-13 DJF # ------------ # # + Added decimal source number and data source filename excluding leading directory to title of spectral plots # # version 1.17 - 2006-06-01 ACS # ------------ # # + Removed filename from label in spectrum plot # # version 1.16 - 2006-04-11 DJF # ------------ # # + Changed use of echeckregion. We cannot use the return code # so the task must write output to a file and we must parse that file. # # version 1.15 - 2006-03-31 DJF # ------------ # # + Modified search for canned rmf. Fail if one is not found. # + Use of TLMIN/TLMAX can leave light curves with all zero 'RATE'. Now check for this and avoid doing anything further with them. # + EXPOSU?? sort was looking for CCD 1 .. 12 for both MOS and PN. Now restricted to 1 .. 7 for MOS. # + Test of BKGDSCRN for flare background screening was testing if keyword was set or not, not if it were T or F (or Y/N ...) # # version 1.14 - 2006-03-15 DJF # ------------ # # + Error check in ML_ID_SRC selection called too early. Tis stopped the test completing. Moved. # # version 1.13 - 2006-03-15 DJF # ------------ # # + Removed attempt to make SSP for PN small window mode. This needs more careful consideration. Without PN SW mode source detection it can be impossible to match SSP selected sources with emldetect source list entries. # # version 1.12 - 2006-03-10 DJF # ------------ # # + fmodhead template creation moved out of if(){} statement scope # so it can be used for additional keyword setting. # # version 1.11 - 2006-03-09 ACS # ------------ # # + Changed writing of source-specific keywords into the header of the spectrum and arf-files # from 'addattribute' to 'fmodhead'; (also corrected keyword comments for addattribute in @kwdCommentList) # # version 1.10 - 2006-03-03 ACS # ------------ # # + Added very crude solution to extract source-level products for PN SW when the source lists # have no PN entries at all (and esources therefore indicates that no PN sources # will be extracted). Now _all_ sources in the SSP source list will be used, # and echeckregion will through out all sources _not_ in the detector window. # # version 1.09 - 2006-03-02 ACS # ------------ # # + Removed sections about TSART and TSTOP keywords (see v 1.07) # since it was done in the wrong place with the wrong values; # moved this part to EPICSources.pm # # version 1.08 - 2006-03-02 ACS # ------------ # # + Add echeckregion after making source extraction region, # will skip the source if it lies outside the detector window # # version 1.07 - 2006-02-21 DJF # ------------ # # + Write TSART and TSTOP keywords from esources output to # lccorr output as TLMIN1 and TLMAX1 ahead of calling evselect. # This is to cause timeseries to be aligned on these values. # # # version 1.06 - 2006-02-20 DJF # ------------ # # + Sources were being selected with an foo() and info() if X # This depends on info() returning true which is not guarenteed. # changed to a more normal if(X) { foo(); info(); } # # version 1.05 - 2006-02-07 DJF # ------------ # # + Sort events for MOS and PN. There are a few ODF where MOS frames are not in order. # This causes lccorr problems so we sort them first. # # version 1.04 - 2006-01-25 DJF # ------------ # # + Fixed unnecessary (and problematic) quoting of ftool parameters. # # version 1.03 - 2005-12-06 DJF # ------------ # + Raise an exception if catASCIIFiles fails. # # # version 1.02 - 2005-11-30 IMS # ------------ # + Added --maketimecolumn=yes and --makeratecolumn=yes to the new evselect calls. # # version 1.01 - 2005-11-30 IMS # ------------ # + Changed single quotes to double in --expression values for the two new evselect calls. # # version 1.0 - 2005-11-29 IMS # ------------ # + I've replaced etimeget by 2 calls of evselect. This is because I find that etimeget does filtering on FLAG and PATTERN without being told to. We need to keep control of this within the pipeline. # + SSPs now not made unless the exposure was flare-screened (requested by ACS 20051129). # + hasFITSKeyword added to the list of used subroutines of ModuleResources. # # version 0.24 - 2005-11-25 IMS # ------------ # + --updateexposure=no added to 2nd evselect call. # # version 0.23 - 2005-11-25 IMS # ------------ # + Existence of spectrum PS file now checked before PStoPDF is invoked. # + I'd forgotten to copy the sorted event list back to $forTimeFilteredEvList. D'oh!! # + The DS9 region product is now not made until after the time series and the spectrum. There's no point in making a DS9 region if there is neither useful spectral or TS data. # # version 0.22 - 2005-11-24 IMS # ------------ # + TS products now labelled as band 8. # # version 0.21 - 2005-11-23 IMS # ------------ # # + Replaced the first evselect call (filters the ml source list) by fselect. This is because evselect failed with error 'Accessing data SRCLIST.CCDNR of type Int32 as int8 is not supported'. # # version 0.20 - 2005-11-23 IMS # ------------ # # + Altered it so that the emldetect source list (suitably filtered) is supplied to region instead of the SSP-filtered list (made within EPICSources from the srcmatch/eposcorr/evalcorr output list). This requires some careful accounting to make sure that we are always talking about the same source. # + Now tests for existence of an EXPOSUnn extension before attempting to sort it (see ChangeLog entry 0.18). # + Added hasFITSExtension to the list of subroutines called from ModuleResources. # # version 0.19 - 2005-11-22 IMS # ------------ # # + Changed the order of the ds9 components in catASCIIFiles so the white circles don't get overdrawn with cyan. # # version 0.18 - 2005-11-21 IMS # ------------ # # + Module now tests that time series have more than 0 non-null RATE values (efftplot etc not called if not). # + EXPOSU extensions of pn event lists are now sorted. # # version 0.17 - 2005-11-16 IMS # ------------ # # + Added a couple more infos. # # version 0.16 - 2005-11-16 IMS # ------------ # # + Module no longer fails if efftplot fails - just the subsequent psToPdf is not called. # # version 0.15 - 2005-11-15 IMS # ------------ # # + Added a print of the rmf file name as read from the temp spectrum, in an attempt to track down why the module can't work out the path name. # + Added carriage returns to the fmodhead template file! # # version 0.14 - 2005-11-14 IMS # ------------ # # + Incorporating slight changed DJF made to intermediate file names. # + Somewhat more care is now taken in juggling the ARF, RMF and bkgSpectrum file vs path names in the source-spectrum SPECTRUM extension header. The point is that xspec needs the path name; but after the xspec call, because pcms path names are meaningless to the user, just the bare file names are substituted, using now a second grppha call. # # version 0.13 - 2005-11-03 IMS # ------------ # # + Label of spectrum pdf plot set to 'source minus background' instead of the default. # + fmodhead template file (introduced in v0.9) was left with just place-holder values! These have now been corrected. # # version 0.12 - 2005-10-31 IMS # ------------ # # + Reversed extraction of bare file name from the path name for the rmf file. I think especget writes this bare name to the header already; don't have to fiddle with it. But check this. # + Changed content string for source ds9 region product to be consistent with FileMap. # # version 0.11 - 2005-10-28 IMS # ------------ # # + Modified the xspec command file a little in order to change the y-axis scale of spectral plots from linear to log (request from ACS). # + Changed the grppha invocation such that the bare file names rather than path names of the arf, rmf and bkg spectrum files are written as keywords to the src spectrum headers. # # version 0.10 - 2005-10-27 IMS # ------------ # # + Now setting the values of lccorr_pcms parameters srcra and srcdec; choice for --srcweightstyle now changed from 'events' to 'deduce'; new parameter --srcdeducestyle='events' now set. # # version 0.09 - 2005-10-14 IMS # ------------ # # + Shortened the content string for the SSP source list. # + Call to lccorr_pcms now uses --srcweightstyle='events' rather than the default 'psf'. It is hoped that this will cure lccorr_pcms's remaining seg fault problems. # + Added numberFITSRows, readFITSKeyword to list of used subroutines of ModuleResources. # + elcplot and efftplot are now run on $correctedRateSet (the final product TS) instead of $bkgSubtractedRateSet, which is an intermediate file. # + The first call to addattribute (the one which writes to the time series) has temporarily been replaced by a call to fmodhead. This is because, whereas elcplot (correctly) expects these keywords to be in the RATE extension header, addattribute is not yet capable of writing them to anywhere but the primary header. (Note that this ought not to stuff up efftplot since this task looks in both places for the kwds.) # + Added new product (ds9 source regions). # + Replaced a few instances of $$srcIdCol[$i] by $id. # # version 0.08 - 2005-09-23 RGW # ------------ # # + reinstate chi-square test in ekstest # # # version 0.07 - 2005-09-23 IMS # ------------ # # + Various kwds get written to the SSP files. However they were all being written to :PRIMARY. These calls to addattribute have now been changed so that the kwds are written to the header of the main binary-table extension (= RATE for light curves, SPECTRUM for spectra and SPECRESP for the arf). # + The RA_CORR and DEC_CORR columns are now being read from the $sspSrcFile. Previously the RA and DEC values were substituted. # + Call to efftplot was commented out (?) - now restored. # # version 0.06 - 2005-09-14 IMS # ------------ # # + Added parameters --tempset and --tempeventset to lccorr_pcms call. # + Failure in lccorr_pcms now does not cause the module to fail. # # version 0.05 - 2005-09-14 IMS # ------------ # # + Added extra parameters --srcregiontab and --bkgregiontab to the call to lccorr_pcms. # # version 0.04 - 2005-09-09 RGW # ------------ # # + fixed arguments to addattribute (integer, not int) # + specify GTI block name to gtialign # + use pre-defined content name for time series & plots # + work around short filename limits in XSPEC # + grppha requires "!" on output file to force overwrite # + fixed parameter names on especget # + not necessary to pass src number as a hex string to newFile # + added a trap to test for lightcurves with zero bins (which can # happen in observations which contain both full-frame and windowed # exposures) # + complete gracefully when no sources have been selected to product creation # + added trap to ignore empty flare GTIs # # Version 0.03 - 2005-09-05 IMS # ------------ # # + Now calling lccorr_pcms, with 8 additional parameters recording the event selections, instead of lccorr. # # version 0.02 - 2005-08-24 RGW # ------------ # # + temporary hacks to get it running within the PCMS # # version 0.01 - 2005-06-09 IMS # ------------ # # + First draft. Anja's script converted to module format. # # Declare list of instruments this module is interested in @instList = qw(epn emos1 emos2); # Number of streams sub numberOfStreams { return numberOfExposures(); } # Rules method sub evaluateRules { # If upstream module has been ignored, so if this one return ignore() if ignored(module => 'ExpChooser' ,instrument => thisInstrument ,stream => thisStream ); return ignore() if ignored( module => 'EPICSources' , instrument => 'all' , stream => 1 ); return start() if complete( module => 'EPICSources' , instrument => 'all' , stream => 1 ); } # Action method sub performAction { info("Module version number: $version"); my $exp_id = exposureID( instrument => thisInstrument , stream => thisStream ); info("Exposure ID: $exp_id"); my $flareScreened = undef; # Make sure there are no empty band images and flare background screening has been applied. my @imgcount; foreach my $band ( 1 .. 5 ) { my $image = findFile( class => 'product' , instrument => thisInstrument , exp_id => $exp_id , band => $band , content => 'EPIC image' , 'format' => 'FITS' ); next unless ($image && fileExists(file =>$image)); if (!defined($flareScreened)) { if ( hasFITSKeyword( file => $image , extension => 'PRIMARY' , keyword => 'BKGDSCRN' )) { $flareScreened = readFITSKeyword( file => $image , extension => 'PRIMARY' , keyword => 'BKGDSCRN' ); $flareScreened = ($flareScreened =~ /$RE{boolean}{true}/i) ? 1 : 0 ; if (defined $flareScreened && !$flareScreened) { info("Exposure was not flare screened"); info("Source-specific products will not be made."); return success(); } } } my @fopt = ( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , format => "ASCII" , band => $band , content => "Filtered image statistics" ); my $stats = findFile( @fopt ); unless ($stats ) { $stats = newFile( @fopt ); } if ( ModuleUtil::isImageEmpty($image, $stats) ) { info("$image is empty"); info("Excluding this exposure from source specific product creation."); return success(); } push @imgcount , $image; } unless (@imgcount) { info("This exposure has no images."); info("Excluding this exposure from source specific product creation."); return success(); } # Filter the ml source list: # my $mlList = findFile( class => 'product' , instrument => 'epic' , content => 'EPIC observation ml source list' ); return success() unless $mlList; my $selExpression; if (thisInstrument eq 'emos1') { $selExpression = "(DET_ML>15 && (ID_BAND==0)&&(ID_INST==2))"; } elsif (thisInstrument eq 'emos2') { $selExpression = "(DET_ML>15 && (ID_BAND==0)&&(ID_INST==3))"; } else # assume 'epn' { $selExpression = "(DET_ML>15 && (ID_BAND==0)&&(ID_INST==1))"; } my $filteredMlList = newFile( class => 'intermediate' , instrument => thisInstrument , content => 'Filtered ML src list' ); doCommand( 'fselect' , infile => $mlList."[SRCLIST]" , outfile => $filteredMlList , expr => $selExpression ) or return exception(); my $mlSrcIdColInOrig = readFITSColumn( file => $filteredMlList , extension => "SRCLIST" , column => "ML_ID_SRC" ); # Find the source-specific-products source list: # my $sspSrcFile = findFile( class => 'intermediate' , instrument => 'epic' , content => 'SSP source list' ); return success() if ( !defined($sspSrcFile) || $sspSrcFile eq '' || !fileExists( file => $sspSrcFile ) ); my $nsrc=numberFITSRows(file => $sspSrcFile, extension => "SRCLIST"); info("Source list contains $nsrc sources"); return success() if $nsrc==0; my $flareGtifile = findFile( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , content => 'flare GTI' ); my $flareGtiExists = 0; if (fileExists(file => $flareGtifile) && numberFITSRows(file => $flareGtifile, extension => "STDGTI")) { $flareGtiExists = 1; } my $forTimeFilteredEvList = findFile( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , content => 'Event list for light curves' ); if ( !defined($forTimeFilteredEvList) || $forTimeFilteredEvList eq '' || !fileExists( file => $forTimeFilteredEvList ) ) { info('No light curve event list. Nothing to do.'); return success(); } ### why is this 'success' but for the spectrum evlist 'exception'?? my $forSpecFilteredEvList = findFile( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , content => 'Event list for spectra' ); unless ( $forSpecFilteredEvList && fileExists( file => $forSpecFilteredEvList ) ) { info('No spectra event list. Nothing to do.'); return success(); } my $flagBitKwdName = 'PN_P_BIT'; if (thisInstrument eq 'emos1') { $flagBitKwdName = 'M1_P_BIT'; } elsif (thisInstrument eq 'emos2') { $flagBitKwdName = 'M2_P_BIT'; } my $flagBit = -1; # default $flagBit = readFITSKeyword( file => $sspSrcFile , extension => "SRCLIST" , keyword => $flagBitKwdName ); return exception() if ($flagBit < 0); my $flagMask = 1 << $flagBit; # Set up the event selection parameter values for lccorr_pcms: my ($patternHi, $eventFlagMask, $rawylo, @numCcds); if (thisInstrument eq 'epn') { $patternHi = 4; $eventFlagMask = '0xfffffef'; $rawylo = 0; @numCcds = ( 1 .. 12 ); } else # assume mos { $patternHi = 12; $eventFlagMask = '0x766ba000'; $rawylo = 0; @numCcds = ( 1 .. 7 ); } # NOTE! The flag masks should be the same as the ones actually used in EPICSources. my @gtiTabList; foreach my $i (@numCcds) { push @gtiTabList, sprintf('%s:STDGTI%02d' , $forTimeFilteredEvList , $i); } my ($timeBinSizeCol, $srcIdCol, $srcIdHexCol, $mlSrcIdColInSSP); my ($srcRaCol, $srcDecCol, $srcCorrRaCol, $srcCorrDecCol); my @selectedSources; if ($sspSrcFile) { $timeBinSizeCol = readFITSColumn( file => $sspSrcFile , extension => "SRCLIST" , column => "TBIN_WIDTH" ); $srcIdCol = readFITSColumn( file => $sspSrcFile , extension => "SRCLIST" , column => "SRC_NUM" ); $srcIdHexCol = readFITSColumn( file => $sspSrcFile , extension => "SRCLIST" , column => "SRC_NUM_HEX" ); $mlSrcIdColInSSP = readFITSColumn( file => $sspSrcFile , extension => "SRCLIST" , column => "ML_ID_SRC" ); $srcRaCol = readFITSColumn( file => $sspSrcFile , extension => "SRCLIST" , column => "RA" ); $srcDecCol = readFITSColumn( file => $sspSrcFile , extension => "SRCLIST" , column => "DEC" ); $srcCorrRaCol = readFITSColumn( file => $sspSrcFile , extension => "SRCLIST" , column => "RA_CORR" ); $srcCorrDecCol = readFITSColumn( file => $sspSrcFile , extension => "SRCLIST" , column => "DEC_CORR" ); my $flagCol = readFITSColumn( file => $sspSrcFile , extension => 'SRCLIST' , column => 'PROCESS' ); for ( my $i=0 ; defined $$flagCol[$i] ; $i++ ) { return exception() if ( ! ( defined($$timeBinSizeCol[$i]) && defined($$srcIdCol[$i]) && defined($$srcIdHexCol[$i]) && defined($$mlSrcIdColInSSP[$i]) )); if (($$flagCol[$i] & $flagMask) > 0) { push( @selectedSources, $i ); info("select source $i"); } } } my %mlId_to_rowNum = (); foreach my $i (@selectedSources) { #info("row check i=$i , \$mlSrcIdColInSSP->[$i] = $mlSrcIdColInSSP->[$i] , j = ( 0 .. $#$mlSrcIdColInOrig ) "); foreach my $j (0 .. $#$mlSrcIdColInOrig) { #info("row check j=$j , \$mlSrcIdColInOrig->[$j] = $mlSrcIdColInOrig->[$j] "); if ($mlSrcIdColInOrig->[$j] == $mlSrcIdColInSSP->[$i]) { $mlId_to_rowNum{ $mlSrcIdColInSSP->[$i] } = $j+1; } } unless ( exists( $mlId_to_rowNum{ $mlSrcIdColInSSP->[$i] } ) && defined( $mlId_to_rowNum{ $mlSrcIdColInSSP->[$i] } ) ) { info("Can't find ML_ID_SRC value $mlSrcIdColInSSP->[$i] in $filteredMlList"); return exception(); } } ## ## END: 2006-03-15 , DJF ## # Sort the EXPOSU extensions of event lists. # This has to be done for both PN and MOS my $tempEventList = newFile( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , content => 'Temporary event list' ); copyFile(source => $forTimeFilteredEvList, destination => $tempEventList); foreach my $i (@numCcds) { my $extname = sprintf('EXPOSU%02d', $i); if (!hasFITSExtension(file => $tempEventList, extension => $extname)) { info("File $tempEventList doesn't contain extension $extname - sort won't be performed."); next; } doCommand( 'fsort' , infile => $tempEventList."[$extname]" , columns => 'TIME' , method => 'insert' ) or return exception(); } copyFile(source => $tempEventList, destination => $forTimeFilteredEvList); # Now make the products: # Record the fact this exposure is used for source specific products in the database setExposureProperty( instrument => thisInstrument , exp_id => $exp_id , name => 'ssp' , value => 0+TRUE ); my $allSrcRegionFile = newFile( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , content => 'all src ds9 region' , format => 'ASCII' ); doCommand( 'slconv' , srclisttab => "${filteredMlList}:SRCLIST" , radiusexpression => 60 , radiusstyle => 'raw' , withlabels => 'no' , colour => 'cyan' , outputstyle => 'ds9' , outfilestyle => 'whole' , outfile => $allSrcRegionFile ) or return exception(); foreach my $i (@selectedSources) { my @kwdNameList = qw(SRCNUM SRC_RA SRC_DEC SRC_CRA SRC_CDEC); my @kwdTypeList = qw(integer real real real real); my @kwdRealValsList = ($$srcRaCol[$i], $$srcDecCol[$i] , $$srcCorrRaCol[$i], $$srcCorrDecCol[$i]); my @kwdWithCommentList = qw(y y y y y); my @kwdCommentList = ("EPIC source list number" , "RA of source", "DEC of source", "Corrected RA of source" , "Corrected DEC of source"); #### units not yet available as a parameter to addAttribute? my $id=$$srcIdCol[$i]; my $src_ra=$$srcRaCol[$i]; my $src_dec=$$srcDecCol[$i]; my $row_in_filteredML=$mlId_to_rowNum{ $mlSrcIdColInSSP->[$i] }; my $srcRegionSet1 = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Source region set 1' ); my $srcRegionSet2 = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Source region set 2' ); my $bkgRegionSet1 = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Background region set 1' ); my $bkgRegionSet2 = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Background region set 2' ); my $tempSrcList = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Temporary source list' ); doCommand( 'region' , eventset => $forTimeFilteredEvList , srclisttab => "${filteredMlList}:SRCLIST" , operationstyle => 'single' , usersrcid => $row_in_filteredML , regionset => $srcRegionSet1 , bkgregionset => $bkgRegionSet1 , expression => '' , radiusstyle => 'userfixed' , fixedradius => 28.0 ## -> variable in the script? , shrinkconfused => 'no' , outunit => 'xy' , tempset => $tempSrcList ) or return exception(); doCommand( 'region' , eventset => $forTimeFilteredEvList , srclisttab => "${filteredMlList}:SRCLIST" , operationstyle => 'single' , usersrcid => $row_in_filteredML , regionset => $srcRegionSet2 , bkgregionset => $bkgRegionSet2 , expression => '' , radiusstyle => 'userfixed' , fixedradius => 60.0 ## -> variable in the script? , shrinkconfused => 'no' , outunit => 'xy' , tempset => $tempSrcList ) or return exception(); # Check whether source lies within the detector window with Richard's new # task 'echeckregion' in the especget package; should replace the subsequent # check on the bins of the light curve (see below) but keep that one in as # a safe guard my $echeckregion_output = newFile(class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , src_num => $id , content => 'echeckregion output' , format => 'ASCII' ); my $echeckregion_status = ModuleUtil::isSourceOnChip( "region($srcRegionSet1,X,Y)" , $forTimeFilteredEvList , $echeckregion_output ); if (!$echeckregion_status) { info("Source lies outside detector window, skipping this source for this exposure."); next; } else { info("Source seems to be within detector window. Continue making products."); } # Make and process the light curve for this source: my $srcRateSet = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Uncorrected source time series' ); my $bgdRateSet = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Uncorrected background time series' ); doCommand( 'evselect' , table => "${forTimeFilteredEvList}:EVENTS" , writedss => 'yes' , updateexposure => 'yes' , expression => "region($srcRegionSet1,X,Y)" , energycolumn => 'PI' , withrateset => 'yes' , rateset => $srcRateSet , timebinsize => $$timeBinSizeCol[$i] , maketimecolumn => 'yes' , makeratecolumn => 'yes' ) or return exception(); doCommand( 'evselect' , table => "${forTimeFilteredEvList}:EVENTS" , writedss => 'yes' , updateexposure => 'yes' , expression => "region($bkgRegionSet2,X,Y)" , energycolumn => 'PI' , withrateset => 'yes' , rateset => $bgdRateSet , timebinsize => $$timeBinSizeCol[$i] , maketimecolumn => 'yes' , makeratecolumn => 'yes' ) or return exception(); # Check whether the lightcurve has any bins, otherwise skip any further # steps. This catches cases where, eg. full-frame imaging overlaps window # mode data (so sources exist in the list where there are no photons # in the windowed exposure) my @rate = readFITSColumn(file => $srcRateSet , extension => 'RATE' , column => 'RATE' ); unless(grep { $_ } @rate) { info("Lightcurve contains no usefull bins, skipping this source for this exposure"); next; } my $bkgSubtractedRateSet = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Exposure-corrected time series' ); my $tempTS = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Temporary time series' ); my $tempEventList = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Temporary event list' ); my $lccorr_succeeded = doCommand( 'lccorr_pcms' , srctsset => $srcRateSet , eventset => $forTimeFilteredEvList , withbkgtsset => 'yes' , bkgtsset => $bgdRateSet , subtractbkg => 'yes' , treatvignet => 'no' , treatfiltertrans => 'no' , treatquantumeff => 'no' , treatcosmicrays => 'yes' , treatdeadtime => 'yes' , treatgti => 'yes' , treatalias => 'no' , applycorrections => 'yes' , outset => $bkgSubtractedRateSet , srcweightstyle => 'deduce' # this choice not supported in 'classical' lccorr. , srcra => $src_ra , srcdec => $src_dec , srcdeducestyle => 'events' # this parameter not supported in 'classical' lccorr. # The following are parameters added to lccorr_pcms to compensate for it not reading the DSS of the rate file $srcRateSet. It is up to YOU to make sure that the parameter values correctly map the rate file event selections. , patternlo => 0 , patternhi => $patternHi , flagmask => $eventFlagMask , pilos => 200 , pihis => 12000 , withgtis => 'yes' , gtitabsnode0 => [@gtiTabList] , rawylo => $rawylo , srcregiontab => "$srcRegionSet1:REGION" , bkgregiontab => "$bkgRegionSet2:REGION" , tempset => $tempTS , tempeventset => $tempEventList ); my $process_ts_further = 0; # default if (!$lccorr_succeeded || !fileExists(file => $bkgSubtractedRateSet)) { info("Lccorr failed for source $id"); # Skip rest of TS processing for this source, go straight to # start making spectral products. } else { # lccorr succeeded, process further: # I don't think this is needed anymore. # doCommand( # 'evselect' # , table => $bkgSubtractedRateSet # , keepfilteroutput => 'yes' # , withfilteredset => 'yes' # , filteredset => $tempTS # , writedss => 'no' # , updateexposure => 'no' #FISH - spong evselect and n columns test below # , expression => '!(isNull(RATE) && RATE > 0)' # ) or return exception(); # Make sure the time series has non-null/zero bins before continuing. my @rate = readFITSColumn(file => $bkgSubtractedRateSet , extension => 'RATE' , column => 'RATE' ); if (grep { $_ } @rate) { $process_ts_further = 1; } else { info("Time series has no rows with non-null values of RATE. Not processed further."); } } if ($process_ts_further) { doCommand( 'addattribute' ,attributename => [qw(SRCNUM SRC_RA SRC_DEC SRC_CRA SRC_CDEC)] ,attributetype => [qw(integer real real real real)] ,attributecomment => [ 'EPIC source list number' ,'RA of source' ,'DEC of source' ,'Corrected RA of source' ,'Corrected DEC of source' ] ,integervalue => $id ,realvalue => [ $$srcRaCol[$i] , $$srcDecCol[$i] , $$srcCorrRaCol[$i] , $$srcCorrDecCol[$i] ] , set => "$bkgSubtractedRateSet:RATE" ) or exception(); my $mergedGtis = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Merged GTIs' ); my @gtiTableList = ("${bkgSubtractedRateSet}:SRC_GTIS" , "${bkgSubtractedRateSet}:BKG_GTIS"); if ($flareGtiExists) {push @gtiTableList, "${flareGtifile}:STDGTI";} doCommand( 'gtimerge' , tables => [@gtiTableList] , withgtitable => 'yes' , gtitable => $mergedGtis , mergemode => 'and' ) or return exception(); my $alignedGtis = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Aligned merged GTIs' ); doCommand( 'gtialign' , style => 'generic' , ingtitable => "$mergedGtis:STDGTI" , outgtitable => "$alignedGtis:STDGTI" , tstable => "${bkgSubtractedRateSet}:RATE" ) or return exception(); my $correctedRateSet = newFile( class => 'product' , instrument => thisInstrument , band => 8 , src_num => $id , exp_id => $exp_id , content => 'EPIC source timeseries' ); doCommand( 'ekstest' , set => $bkgSubtractedRateSet , gtis => 'yes' , gtiset => $alignedGtis , chisquaretest => 'yes' , screen => 'yes' , newoutset => 'yes' , outset => $correctedRateSet ) or return exception(); my $psFile = newFile( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , src_num => $id , content => 'Source time series plot' , format => 'PS' ); my $pdfFile = newFile( class => 'product' , instrument => thisInstrument , band => 8 , exp_id => $exp_id , src_num => $id # newFile takes care of hexification , content => 'EPIC source timeseries plot' , format => 'PDF' ); doCommand( 'elcplot' # , set => $bkgSubtractedRateSet , set => $correctedRateSet , usegtiset => 'yes' , useflareset => 'no' , gtiset => $alignedGtis , offset => 'yes' , offsetstyle => 'auto' , forcestart => 'yes' , binsize => 1 , plotdevice => '/vcps' , plotfile => $psFile ) or return exception(); PStoPDF( source => $psFile , destination => $pdfFile ) or return exception(); my $fftPsFile = newFile( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , src_num => $id , content => 'Time series FFT plot' , format => 'PS' ); my $fftPdfFile = newFile( class => 'product' , instrument => thisInstrument , band => 8 , exp_id => $exp_id , src_num => $id # newFile takes care of hexification , content => 'EPIC source FFT plot' , format => 'PDF' ); # doCommand( my $efftplot_succeeded = doCommand( 'efftplot' , infile => $correctedRateSet , xscale => 'log' , yscale => 'lin' , gtis => 'yes' , gtiset => $alignedGtis , plotdevice => '/vcps' , outfile => $fftPsFile ); # or return exception(); if (!$efftplot_succeeded || !fileExists(file => $fftPsFile)) { info("efftplot failed for source $id"); } else { PStoPDF( source => $fftPsFile , destination => $fftPdfFile ) or return exception(); } } # end 'if lccorr worked' sequelae # Now create spectral products: my $srcSelExpression = "(region($srcRegionSet1,X,Y))"; my $bkgSelExpression = "(region($bkgRegionSet2,X,Y))"; if ($flareGtiExists) { $srcSelExpression .= " && (GTI(${flareGtifile}:STDGTI,TIME))"; $bkgSelExpression .= " && (GTI(${flareGtifile}:STDGTI,TIME))"; } my $tempSrcSpectrum = newFile( class => 'intermediate' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'Source spectrum' ); my $bkgSpectrum = newFile( class => 'product' , instrument => thisInstrument , src_num => $id # newFile takes care of hexification , exp_id => $exp_id , content => 'EPIC source background spectrum' ); my $arf = newFile( class => 'product' , instrument => thisInstrument , src_num => $id # newFile takes care of hexification , exp_id => $exp_id , content => 'EPIC ancillary response function' ); doCommand( 'especget' , table => $forSpecFilteredEvList , srcexp => $srcSelExpression , backexp => $bkgSelExpression , withfilestem => 'no' , srcspecset => $tempSrcSpectrum , bckspecset => $bkgSpectrum , srcarfset => $arf , witharfset => 'yes' , withrmfset => 'yes' # but leave --rmfset at default (yeesh, # what an interface!) - this makes the # task look for a canned rmf. , useodfatt => 'no' , withbadpixcorr => 'yes' , extendedsource => 'no' ) or return exception(); # especget in this configuration doesn't either create or # need an RMF file: all it does is decide which of the canned # RMF matrices at the soc (sort of a pseudo-ccf) is the # correct one to use for this source, filemode etc. Especget # writes the name of this file (just the file name, no path # or url information) to the keyword RESPFILE in the source # spectrum dataset header. Apparently this keyword is where # xspec expects to find the RMF file name. # # This whole section is a bit ludicrous. All that is required # is to make a ps plot of the spectrum. There was a sas task, # especplot, designed to do this. However there seemed to be # eternal problems with especplot. For political reasons it # was deemed inadvisable for an alternate developer to code # this task. So it was decided to use xspec to make this plot. # However xspec requires an RMF. Hence it is necessary to keep # a cache of ALL the canned rmfs - just to make a plot!! # # In order for the whole Heath Robinson contraption to work, # the simple file name in the RESPFILE keyword must be # overwritten by a string of the form private_pcms_rmf_path/ # filename. The the simple file name must be put back into # RESPFILE, for the benefit of a later interactive user of # the spectrum. Eesh what a horrid system. # Enter source-specific keywords into spectral files: doCommand( 'addattribute' ,attributename => [qw(SRCNUM SRC_RA SRC_DEC SRC_CRA SRC_CDEC)] ,attributetype => [qw(integer real real real real)] ,attributecomment => [ 'EPIC source list number' ,'RA of source' ,'DEC of source' ,'Corrected RA of source' ,'Corrected DEC of source' ] ,integervalue => $id ,realvalue => [ $$srcRaCol[$i] , $$srcDecCol[$i] , $$srcCorrRaCol[$i] , $$srcCorrDecCol[$i] ] , set => ["$tempSrcSpectrum:SPECTRUM" , "$bkgSpectrum:SPECTRUM", "$arf:SPECRESP" ] ) or exception(); my $srcSpectrum = newFile( class => 'product' , instrument => thisInstrument , src_num => $id , exp_id => $exp_id , content => 'EPIC source spectrum' ) or exception(); # For each of the three files ARF, RMF and bkgSpectrum, we need to know both the path name (to pass to xspec, so it can find the files) and the bare file name (to write to the output spectrum). The ARF and bkgSpectrum path names we know already, they are contained respectively in the variables $arf and $bkgSpectrum. These have to be cut up to yield the bare file name. For the RMF on the other hand we can read the bare file name from the SPECTRUM header of $tempSrcSpectrum, where especget writes it. There is a special pcms routine findRmfFile to transform this into a full path name. my $cannedRmfFileName = readFITSKeyword( file => $tempSrcSpectrum , extension => 'SPECTRUM' , keyword => 'RESPFILE' ); info("Canned RMF is $cannedRmfFileName"); # Find the local path name $cannedRmfPathName: my $cannedRmfPathName = findRmfFile( file => $cannedRmfFileName ); unless ( $cannedRmfPathName ) { exception('MISSINGRMF','Unable to locate requested rmf : ', $cannedRmfFileName ); } info("Absolute path to canned RMF is $cannedRmfPathName"); # Before running xspec we need to make sure make sure that # $cannedRmfPathName exists and is readable. Otherwise xspec # will end in an endless loop. - DJHF This is done by findRmfFile() # ie. it returns undef when it cannot find a FILE. # Cut off the directory names from the ARF and bkgSpectrum files: my @chunks; @chunks = split('\/', $arf); my $arfFileName = $chunks[$#chunks]; @chunks = split('\/', $bkgSpectrum); my $bkgSpectrumFileName = $chunks[$#chunks]; my $grpphaCommExpr; if (thisInstrument eq 'epn') { $grpphaCommExpr = "bad 0-63&bad 1575-1645&group min " ."20&chkey RESPFILE $cannedRmfPathName"."&exit"; } elsif (thisInstrument eq 'emos1' || thisInstrument eq 'emos2') { $grpphaCommExpr = "bad 0-23&group min " ."20&chkey RESPFILE $cannedRmfPathName"."&exit"; } # Hope I've at last got this right! IMS 20051114 doCommand( 'grppha' , infile => $tempSrcSpectrum , outfile => "\!$srcSpectrum" , chatter => 0 , comm => $grpphaCommExpr ) or return exception(); # xspecIIIIEEEESSS.ps my $xspecname = sprintf( '%s%s%s%03X' ,'xspec' , thisInstrument , $exp_id , $id); my $spectrumPsFile = newFile( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , src_num => $id , content => 'Spectrum plot' , format => 'PS' , name => "$xspecname.ps" # XSPEC doesn't like long filenames ); my $spectrumPdfFile = newFile( class => 'product' , instrument => thisInstrument , exp_id => $exp_id , src_num => $id # newFile takes care of hexification , content => 'EPIC source spectrum plot' , format => 'PDF' ); my $xpecCommandFile = newFile( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , src_num => $id , content => 'XSPEC command file' , format => 'ASCII' ); # Data source name without leading directories for use in spectrum title. my $fname = (split /\//,$srcSpectrum)[-1]; $fname =~ s/\.FIT$|\.FTZ$//g; my $title = "Source $id minus background ($fname)\n"; open(XSPLIS, ">$xpecCommandFile"); print XSPLIS "data $srcSpectrum\n"; print XSPLIS "ignore **-0.2 12.0-**\n"; print XSPLIS "ignore bad\n"; print XSPLIS "setplot en\n"; print XSPLIS "cpd /NULL\n"; print XSPLIS "setplot command log y on\n"; print XSPLIS "setplot command r y\n"; # print XSPLIS "setplot command ld\n"; print XSPLIS "setplot command r x 0.2 12\n"; print XSPLIS "setplot command lwidth 6\n"; print XSPLIS "setplot command la FILE \n"; print XSPLIS "setplot command la T $title \n"; print XSPLIS "setplot command hardcopy $spectrumPsFile/ps\n"; print XSPLIS "plot\n"; print XSPLIS "exit\n"; print XSPLIS "Y\n"; close(XSPLIS); doCommand( "xspec - '$xpecCommandFile'") or return exception() ; # If the spectrum file SPECTRUM extension has an ANCRFILe keyword # rewrite the FIT to FTZ to maintain post-compression usability. my $ancrfile = readFITSKeyword( file => $srcSpectrum , extension => 'SPECTRUM' , keyword => 'ANCRFILE' ); if ( $ancrfile && $ancrfile =~ /\.FIT/ ) { # Make the keyword refer to a compressed file as this is what the user will # ultimately have. $ancrfile =~ s/\.FIT/\.FTZ/g; doCommand('addattribute' ,set => "$srcSpectrum:SPECTRUM" ,attributename => 'ANCRFILE' ,attributetype => 'string' ,stringvalue => $ancrfile ,attributecomment => "\"\\\"ancillary response\\\"\"" ) or exception(); } # Now rewrite the $cannedRmfFileName back into the RESPFILE keyword!!!! # doCommand( # 'addattribute' # , set => "${srcSpectrum}:PRIMARY" ##### whoops - addattribute is not yet able to write kwds to an extension header! 'PRIMARY' must be wrong in any case. But it seems not to matter here... # , attributename => 'RESPFILE' # , attributetype => 'string' # , stringvalue => $cannedRmfFileName # ) # or return exception(); # INSERTED BY DLG 26.01.09 --------------------- # Manual modification of ANCRFILE and FITSFILE headers # to reflect compressed FITS filenames $arfFileName =~ s/\.FIT/\.FTZ/g; $bkgSpectrumFileName =~ s/\.FIT/\.FTZ/g; $grpphaCommExpr = "chkey ANCRFILE $arfFileName"."&chkey RESPFILE " ."$cannedRmfFileName"."&chkey BACKFILE $bkgSpectrumFileName"."&exit"; # INSERTED BY DLG 26.01.09 --------------------- doCommand( 'grppha' , infile => $srcSpectrum , outfile => "\!$srcSpectrum" , chatter => 0 , comm => $grpphaCommExpr ) or return exception(); if (fileExists(file => $spectrumPsFile)) { PStoPDF( source => $spectrumPsFile , destination => $spectrumPdfFile ) or return exception(); } else { info("Can't find $spectrumPsFile - so can't make PDF."); } if (!(fileExists(file => $spectrumPsFile) || $process_ts_further)) { info("There is neither a PS spectrum product nor a time series with any"); info("useful data; thus DS9-region product for this source won't be made."); next; } # Make the DS9 region product file: my $srcRegionFile = newFile( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , src_num => $id , content => 'src ds9 region' , format => 'ASCII' ); my $bkgRegionFile1 = newFile( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , src_num => $id , content => 'bkg ds9 region 1' , format => 'ASCII' ); my $bkgRegionFile2 = newFile( class => 'intermediate' , instrument => thisInstrument , exp_id => $exp_id , src_num => $id , content => 'bkg ds9 region 2' , format => 'ASCII' ); doCommand( 'slconv' , srclisttab => "${sspSrcFile}:SRCLIST" , expression => "SRC_NUM==$id" , radiusexpression => 28 , radiusstyle => 'raw' , withlabels => 'yes' , labelexpression => 'SRC_NUM' , colour => 'green' , outputstyle => 'ds9' , outfilestyle => 'whole' , outfile => $srcRegionFile ) or return exception(); doCommand( 'slconv' , srclisttab => "${sspSrcFile}:SRCLIST" , expression => "SRC_NUM==$id" , radiusexpression => 60 , radiusstyle => 'raw' , withlabels => 'yes' , labelexpression => 'SRC_NUM' , colour => 'white' , outputstyle => 'ds9' , outfilestyle => 'whole' , outfile => $bkgRegionFile1 ) or return exception(); doCommand( 'slconv' , srclisttab => "${sspSrcFile}:SRCLIST" , expression => "SRC_NUM==$id" , radiusexpression => 180 , radiusstyle => 'raw' , withlabels => 'yes' , labelexpression => 'SRC_NUM' , colour => 'white' , outputstyle => 'ds9' , outfilestyle => 'whole' , outfile => $bkgRegionFile2 ) or return exception(); my $ds9RegionFile = newFile( class => 'product' , instrument => thisInstrument , exp_id => $exp_id , src_num => $id , content => 'EPIC source ds9 region' , format => 'ASCII' ); catASCIIFiles( infiles => [$srcRegionFile, $bkgRegionFile1, $bkgRegionFile2, $allSrcRegionFile] , outfile => $ds9RegionFile ) or return exception(); } # end loop over sources return success(); } 1;