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(:longnames :constants);
use FITSException;
use Background;
use Data::Dumper;
$name = __PACKAGE__;
$VERSION = '1.55';
$version = $VERSION;
$author = 'Anja Schroeder,Ian Stewart,Duncan John Fyfe,Grant Denkinson,Duncan Law-Green,Ed Chapin';
$date = '2013-08-23';
@instList = qw(epn emos1 emos2);
sub numberOfStreams
{
return numberOfExposures();
}
sub evaluateRules
{
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
);
}
sub performAction
{
info("Module version number: $version");
my $exp_id = exposureID( instrument => thisInstrument
, stream => thisStream
);
info("Exposure ID: $exp_id");
my $flareScreened = undef;
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();
}
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
{
$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"
);
my $sspSrcInput = findFile(
class => 'intermediate'
, instrument => 'epic'
, content => 'SSP source list'
);
return success()
if ( !defined($sspSrcInput) || $sspSrcInput eq ''
|| !fileExists( file => $sspSrcInput )
);
my $nsrc=numberFITSRows(file => $sspSrcInput,
extension => "SRCLIST");
info("SSP Source list contains $nsrc sources");
return success() if $nsrc==0;
my $cam;
if (thisInstrument eq 'emos1') {
$cam = 'M1';
} elsif (thisInstrument eq 'emos2') {
$cam = 'M2';
} else {
$cam = 'PN';
};
my $cts_col = $cam."_CTS";
my $offax_col = $cam."_OFFAX";
info("Filtering SSP list $sspSrcInput for $cts_col > 0");
my $sspSrcFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Filtered SSP source list'
, format => 'FITS'
) or return exception();
doCommand(
'fselect'
, infile => $sspSrcInput."[SRCLIST]"
, outfile => $sspSrcFile
, expr => "$cts_col > 0"
);
my $nsrc=numberFITSRows(file => $sspSrcFile,
extension => "SRCLIST");
info("Filtered Source list contains $nsrc sources");
return success() if $nsrc==0;
my $standardFlareGTI = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'flare GTI'
);
my $optFlareGTI = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Optimised flare GTI'
);
my $alignedGtis;
my $flareGtifile;
if (fileExists(file => $optFlareGTI)) {
$flareGtifile = $optFlareGTI;
} else {
$flareGtifile = $standardFlareGTI;
}
my $flareGtiExists = 0;
if (fileExists(file => $flareGtifile) &&
numberFITSRows(file => $flareGtifile,
extension => "STDGTI")) {
$flareGtiExists = 1;
}
my $forTimeFilteredEvList = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Event list for light curves'
);
if ( !defined($forTimeFilteredEvList)
|| $forTimeFilteredEvList eq ''
|| !fileExists( file => $forTimeFilteredEvList )
) {
info('No light curve event list. Nothing to do.');
return success();
}
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;
$flagBit = readFITSKeyword(
file => $sspSrcFile
, extension => "SRCLIST"
, keyword => $flagBitKwdName
);
return exception() if ($flagBit < 0);
my $flagMask = 1 << $flagBit;
my ($patternHi, $eventFlagMask, $rawylo, @numCcds);
if (thisInstrument eq 'epn')
{
$patternHi = 4;
$eventFlagMask = '0xfffffef';
$rawylo = 0;
@numCcds = ( 1 .. 12 );
}
else
{
$patternHi = 12;
$eventFlagMask = '0x766ba000';
$rawylo = 0;
@numCcds = ( 1 .. 7 );
}
my @gtiTabList;
foreach my $i (@numCcds)
{
push @gtiTabList, sprintf('%s:STDGTI%02d' , $forTimeFilteredEvList , $i);
}
my ($timeBinSizeCol, $srcIdCol, $srcIdHexCol, $mlSrcIdColInSSP);
my ($srcRaCol, $srcDecCol, $srcCorrRaCol, $srcCorrDecCol, $srcCtsCol, @srcOffaxCol);
my ($flagCol, $epFlagCol, $pnFlagCol, $m1FlagCol, $m2FlagCol);
my @selectedSources;
if ($sspSrcFile)
{
$timeBinSizeCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "TBIN_WIDTH"
);
$srcIdCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "SRC_NUM"
);
$srcIdHexCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "SRC_NUM_HEX"
);
$mlSrcIdColInSSP = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "ML_ID_SRC"
);
$srcRaCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "RA"
);
$srcDecCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "DEC"
);
$srcCorrRaCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "RA_CORR"
);
$srcCorrDecCol = readFITSColumn(
file => $sspSrcFile
, extension => "SRCLIST"
, column => "DEC_CORR"
);
$flagCol = readFITSColumn(
file => $sspSrcFile
, extension => 'SRCLIST'
, column => 'PROCESS'
);
$epFlagCol = readFITSColumn(
file => $sspSrcFile
, extension => 'SRCLIST'
, column => 'EP_FLAG'
);
$pnFlagCol = readFITSColumn(
file => $sspSrcFile
, extension => 'SRCLIST'
,column => 'PN_FLAG'
);
$m1FlagCol = readFITSColumn(
file => $sspSrcFile
, extension => 'SRCLIST'
, column => 'M1_FLAG'
);
$m2FlagCol = readFITSColumn(
file => $sspSrcFile
, extension => 'SRCLIST'
, column => 'M2_FLAG'
);
$srcCtsCol = readFITSColumn(
file => $sspSrcFile
, extension => 'SRCLIST'
, column => "$cts_col"
);
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)
{
foreach my $j (0 .. $#$mlSrcIdColInOrig)
{
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();
}
}
my ($mlIdCol, $mlRaCol, $mlDecCol, $mlCtsCol, $mlOffaxCol);
my $epicImage = findFile(
class => "product"
, instrument => thisInstrument
, exp_id => $exp_id
, band => 8
, content => 'EPIC image'
, format => 'FITS'
) or return exception();
my $summaryList = findFile(
class => "product"
, instrument => 'epic'
, content => 'EPIC summary source list'
, format => 'FITS'
) or return exception();
my $filteredSrcList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Column filtered source list'
, format => 'FITS'
);
doCommand(
'fselect'
, infile => $summaryList
, outfile => $filteredSrcList
, expr => "$cts_col > 0"
) or return exception();
my $num_sources = numberFITSRows(
file => $filteredSrcList
, extension => 'SRCLIST'
);
info("$filteredSrcList: found $num_sources with ${cam}_CTS > 0");
my $num_ssp = numberFITSRows(
file => $sspSrcFile
, extension => 'SRCLIST'
);
my $subSrcList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Column subset source list'
, format => 'FITS'
) or return exception();
my $columns = "SRC_NUM,$cts_col,RA,DEC,$offax_col";
doCommand(
'ftabcopy'
, infile => $filteredSrcList
, outfile => $subSrcList
, columns => $columns
, rows => "-"
) or return exception();
$mlIdCol = readFITSColumn(
file => $subSrcList
, extension => "SRCLIST"
, column => 'SRC_NUM'
);
$mlRaCol = readFITSColumn(
file => $subSrcList
, extension => "SRCLIST"
, column => "RA"
);
$mlDecCol = readFITSColumn(
file => $subSrcList
, extension => "SRCLIST"
, column => "DEC"
);
$mlCtsCol = readFITSColumn(
file => $subSrcList
, extension => "SRCLIST"
, column => "$cts_col"
);
$mlOffaxCol = readFITSColumn(
file => $subSrcList
, extension => "SRCLIST"
, column => "$offax_col"
);
my ($i, $j);
my (@allX_raw, @allY_raw, @all_CCD);
my $sccd;
for ($i = 0; $i < $num_sources; $i++) {
my $sc = $i+1;
my ($xcoord, $ycoord, $sccd) = coordConv($epicImage, $$mlRaCol[$i], $$mlDecCol[$i]);
info("coordConv returned - source $sc of $num_sources: x: $xcoord y: $ycoord CCD: $sccd");
push @allX_raw, int($xcoord);
push @allY_raw, int($ycoord);
push @all_CCD, int($sccd);
};
my %offax_hash;
@offax_hash{@$mlIdCol} = @$mlOffaxCol;
my @srcCCDCol;
for ( $j = 0; $j < $num_ssp; $j++ ) {
$srcOffaxCol[$j] = $offax_hash{$$srcIdCol[$j]};
$srcCCDCol[$j] = $all_CCD[$$srcIdCol[$j]];
};
my $data;
$$data{X} = [@allX_raw];
$$data{Y} = [@allY_raw];
$$data{-column}{X} = { 'width' => 4,
'colnum' => 6,
'ttype' => 'X',
'repeat' => 1,
'typecode' => 42
};
$$data{-column}{Y} = { 'width' => 4,
'colnum' => 7,
'ttype' => 'Y',
'repeat' => 1,
'typecode' => 42
};
my $colname = $$data{colname};
push (@$colname, 'X', 'Y');
$$data{colname} = [ @$colname ];
writeFITSColumn(
file => $subSrcList,
extension => 'SRCLIST',
column => $data
);
info("Generated X and Y coordinate columns in $subSrcList");
my $lookupFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'Counts radius lookup table'
, format => 'ASCII'
) or return exception();
my $lookupInput = $ENV{PCMS_ROOT} . '/pipeline/etc/LOOKUP.FITS';
copyFile(
source => $lookupInput
, destination => $lookupFile
) or return exception();
info("Copied mask radius lookup table in place");
my $lookupData = readFITSTable(
file => $lookupFile
, extension => 'LOOKUP'
, colname => [qw(CTS_LOW CTS_HIGH RADIUS OFFA_LOW OFFA_HIGH)]
) or return exception();
my $num_lookup = numberFITSRows(
file => $lookupFile
, extension => 'LOOKUP'
);
my @lookup;
for ($i = 0; $i < $num_lookup; $i++) {
push(@lookup, {
'cts_low' => $lookupData->{CTS_LOW}[$i],
'cts_high' => $lookupData->{CTS_HIGH}[$i],
'radius' => $lookupData->{RADIUS}[$i],
'offa_low' => $lookupData->{OFFA_LOW}[$i],
'offa_high' => $lookupData->{OFFA_HIGH}[$i]
});
}
info("Read lookup table $lookupFile...");
my (@allR_raw, @allR_ima);
for ($i = 0; $i < $num_sources; $i++) {
my $cts = $$mlCtsCol[$i];
my $offax = $$mlOffaxCol[$i];
my @result = grep { $_->{cts_low} <= $cts && $_->{cts_high} > $cts &&
$_->{offa_low} <= $offax && $_->{offa_high} > $offax } @lookup;
my $rad = ($result[0]->{'radius'});
exception("Radius undefined for CTS $cts OFFAX $offax") unless defined($rad);
info("Source $i: CTS $cts OFFAX $offax MASK $rad");
$allR_raw[$i] = $rad*20;
$allR_ima[$i] = $rad/4;
}
my $refx1 = readFITSKeyword(
file => $epicImage
, extension => 'PRIMARY'
, keyword => 'REFXCRPX'
) or return exception();
my $refy1 = readFITSKeyword(
file => $epicImage
, extension => 'PRIMARY'
, keyword => 'REFYCRPX'
) or return exception();
my $crpix1 = readFITSKeyword(
file => $epicImage
, extension => 'PRIMARY'
, keyword => 'CRPIX1'
) or return exception();
my $crpix2 = readFITSKeyword(
file => $epicImage
, extension => 'PRIMARY'
, keyword => 'CRPIX2'
) or return exception();
my (@allX_ima, @allY_ima);
for ($i = 0; $i < $num_sources; $i++) {
$allX_ima[$i] = (($allX_raw[$i] - $refx1)/80.0) + $crpix1;
$allY_ima[$i] = (($allY_raw[$i] - $refy1)/80.0) + $crpix2;
}
$$data{RAD_RAW} = [@allR_raw];
$$data{RAD_IMA} = [@allR_ima];
$$data{X_IMA} = [@allX_ima];
$$data{Y_IMA} = [@allY_ima];
$$data{-column}{RAD_RAW} = { 'width' => 4,
'ttype' => 'RAD_RAW',
'repeat' => 1,
'typecode' => 42
};
$$data{-column}{RAD_IMA} = { 'width' => 4,
'ttype' => 'RAD_IMA',
'repeat' => 1,
'typecode' => 42
};
$$data{-column}{X_IMA} = { 'width' => 4,
'ttype' => 'X_IMA',
'repeat' => 1,
'typecode' => 42
};
$$data{-column}{Y_IMA} = { 'width' => 4,
'ttype' => 'Y_IMA',
'repeat' => 1,
'typecode' => 42
};
my $colname = $$data{colname};
pop @$colname;
pop @$colname;
push (@$colname, 'RAD_RAW', 'RAD_IMA', 'X_IMA', 'Y_IMA');
$$data{colname} = [ @$colname ];
writeFITSColumn(
file => $subSrcList,
extension => 'SRCLIST',
column => $data
);
info("Added RAD_RAW, RAD_IMA, X_IMA, Y_IMA columns to $subSrcList");
my $extrRegion = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'extraction region file'
, format => 'FITS'
) or return exception();
my $extrTemp = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'extraction temp file'
, format => 'FITS'
) or return exception();
doCommand(
'fcopy'
, infile => $subSrcList."[1][col X=X_IMA;Y=Y_IMA]"
, outfile => $extrTemp
) or return exception();
doCommand(
'fmerge'
, infiles => $extrTemp
, outfile => $extrRegion
, columns => "X,Y"
, mextname => 'REGION'
) or return exception();
my @SHAPEcol = ("CIRCLE") x $num_sources;
my $data2;
$$data2{SHAPE} = [@SHAPEcol];
$$data2{R} = [@allR_ima];
$$data2{-column}{SHAPE} = { 'width' => 10,
'colnum' => 1,
'ttype' => 'SHAPE',
'repeat' => 10,
'typecode' => 16
};
$$data2{-column}{R} = { 'width' => 4,
'ttype' => 'R',
'repeat' => 1,
'typecode' => 42
};
my $colname2;
push (@$colname2, 'R','SHAPE');
$$data2{colname} = [ @$colname2 ];
writeFITSColumn(
file => $extrRegion,
extension => 'REGION',
column => $data2
);
info("Generated region file $extrRegion");
my $extrDetMask = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'extraction detmask'
, format => 'FITS'
) or return exception();
my $expImage = findFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'EPIC exposure map'
, format => 'FITS'
) or return exception();
doCommand(
'emask'
, expimageset => $expImage
, detmaskset => $extrDetMask
, withregionset => 'yes'
, regionset => $extrRegion
) or return exception();
info("Generated extraction detmask $extrDetMask");
my $pileupMask = findFile(
class => 'intermediate'
, instrument => thisInstrument
, content => 'Pileup mask'
, band => 8
, format => 'FITS'
);
if (fileExists( file => $pileupMask ) ) {
info("Combine extraction detmask with RGA/OOT mask for source extraction");
doCommand(
'farith'
, infil1 => "${extrDetMask}[1]"
, infil2 => $pileupMask
, outfil => $extrDetMask
, clobber => 'yes'
, copyprime => 'yes'
, ops => 'MUL'
) or return exception();
}
my (@X_raw, @Y_raw);
my (@X_ima, @Y_ima);
my (@sCCD);
for ($j = 0; $j < $num_ssp; $j++) {
my $sc = $j+1;
my ($xcoord, $ycoord, $sccd) = coordConv($epicImage, $$srcCorrRaCol[$j], $$srcCorrDecCol[$j]);
$X_raw[$j] = int($xcoord);
$Y_raw[$j] = int($ycoord);
$X_ima[$j] = (($X_raw[$j] - $refx1)/80.0) + $crpix1;
$Y_ima[$j] = (($Y_raw[$j] - $refy1)/80.0) + $crpix2;
$sCCD[$j] = int($sccd);
info("DEBUG: Source $sc of $num_ssp SSPs: X $X_raw[$j] Y: $Y_raw[$j] CCD: $sCCD[$j]");
}
my (@R_raw, @R_ima);
for ($j = 0; $j < $num_ssp; $j++ ) {
my $cts = $$srcCtsCol[$j];
my $offax = $srcOffaxCol[$j];
my @result = grep { $_->{cts_low} < $cts && $_->{cts_high} > $cts &&
$_->{offa_low} < $offax && $_->{offa_high} > $offax } @lookup;
my $rad = ($result[0]->{'radius'});
exception("Radius undefined for CTS $cts OFFAX $offax") unless defined($rad);
$R_raw[$j] = $rad * 20.0;
$R_ima[$j] = $rad / 4.0;
}
info("==== Creating background regions (step 1)");
my $MAX_M = 150;
my @RB_ima = @allR_ima;
my (@XB_raw, @YB_raw, @RB_raw);
my (@XB_ima, @YB_ima);
my $good_area_limit = 0.70;
my $illum = 0.0;
my $max_illum = 0.0;
my $submode = readFITSKeyword(
file => $forSpecFilteredEvList
, extension => 'EVENTS'
, keyword => 'SUBMODE'
) or return exception();
info("DEBUG: Observation submode is $submode");
my $smallWindowMode = 0;
if (thisInstrument =~ /emos/i && $submode =~/PrimePartialW/i) {
$smallWindowMode = 1;
info("Small Window Mode $cam:$submode detected, changing iteration parameters");
}
for ($j = 0; $j < $num_ssp; $j++) {
info("--------------------------------------------------------------");
my $id = $$srcIdCol[$j];
my $sc = $j+1;
my $camcts = $$srcCtsCol[$j];
info("**** Source $id ($sc of $num_ssp): In-camera counts $camcts, CCD: $sCCD[$j]");
info("**** Source $id ($sc of $num_ssp): Step 1 background");
my $m = ($smallWindowMode && $sCCD[$j] == 1) ? 45 : 0;
while ($m < $MAX_M) {
$m = $m + 1;
info("[1] Source ID: $id, iteration $m: Running select_bkg(".$X_ima[$j].",".$Y_ima[$j].",".$R_ima[$j].")");
my $result = &select_bkg($X_ima[$j], $Y_ima[$j], $R_ima[$j], $m, $extrDetMask, $epicImage, $id);
$illum = $$result->{illum};
my $max_index = $$result->{index};
if ($illum > $good_area_limit) {
info("Source $id: CONVERGED @ radial step: $m azimuthal step: $max_index illumination: $illum");
$XB_ima[$j] = $$result->{xback};
$YB_ima[$j] = $$result->{yback};
$RB_ima[$j] = $$result->{radback};
info("Source: $id: background solution XB_ima, YB_ima, RB_ima: ".$XB_ima[$j] . "," . $YB_ima[$j] . "," . $RB_ima[$j]);
last;
} elsif ($illum > $max_illum) {
$XB_ima[$j] = $$result->{xback};
$YB_ima[$j] = $$result->{yback};
$RB_ima[$j] = $$result->{radback};
$max_illum = $illum;
}
}
$XB_raw[$j] = ($XB_ima[$j] - $crpix1) * 80.0 + $refx1;
$YB_raw[$j] = ($YB_ima[$j] - $crpix2) * 80.0 + $refy1;
$RB_raw[$j] = $RB_ima[$j] * 80.0;
$max_illum = 0;
}
info("==== Step 1 iteration complete!");
info("==== Step 2 eregionanalyse optimisation");
my (@X1_raw, @Y1_raw, @R1_raw);
my (@X1_ima, @Y1_ima, @R1_ima);
my $radtmp;
my @source_radius = ();
for ( $j = 0; $j < $num_ssp; $j++ ) {
info("------------------------------------------------------------");
my $id = $$srcIdCol[$j];
my $sc = $j+1;
info("++++ Source ID $id ($sc of $num_ssp): Starting eregionanalyse");
my $cmd = 'eregionanalyse imageset='.$epicImage.' srcexp="(X,Y) IN circle(' .
"$X_raw[$j], $Y_raw[$j], $R_raw[$j])" . '" backexp="(X,Y) IN circle(' .
"$XB_raw[$j], $YB_raw[$j], $RB_raw[$j])" . '"';
info("CMD: $cmd");
my $optoutput = `$cmd`;
info("OPTOUTPUT: $optoutput");
my ($counts) = ($optoutput =~ /counts\sin\ssource\sregion:\s+(\d+)/);
info("Found $counts counts in region");
if ($counts == 0) {
info("Ignoring optimisation...");
$X1_raw[$j] = $X_raw[$j];
$Y1_raw[$j] = $Y_raw[$j];
$R1_raw[$j] = $R_raw[$j];
} else {
($X1_raw[$j]) = ($optoutput =~ /xcentroid:\s+(\d+)/);
($Y1_raw[$j]) = ($optoutput =~ /ycentroid:\s+(\d+)/);
($R1_raw[$j]) = ($optoutput =~ /optradius:\s+\d+\s+arcsecs\s+(\d+)/);
}
if ($R1_raw[$j] < 240.0) {$R1_raw[$j] = 240.0};
if ($R1_raw[$j] > 800.0) {$R1_raw[$j] = 800.0};
$X1_ima[$j] = (($X1_raw[$j] - $refx1)/80.0) + $crpix1;
$Y1_ima[$j] = (($Y1_raw[$j] - $refy1)/80.0) + $crpix2;
$R1_ima[$j] = $R1_raw[$j] / 80.0;
$source_radius[$id] = $R1_raw[$j];
}
info("*** DEBUG DUMP");
for ($j=0; $j < $num_ssp; $j++) {
info("Source $j: raw (".$X1_raw[$j].",".$Y1_raw[$j].",".$R1_raw[$j].") ima (".$X1_ima[$j].",".$Y1_ima[$j].",".$R1_ima[$j].")");
}
info("*** DEBUG DUMP");
info("==== Step 3: Creating background");
my $MAX_M = 150;
my (@XB1_ima, @YB1_ima, @RB1_ima);
my (@XB1_raw, @YB1_raw, @RB1_raw);
for ( $j = 0; $j < $num_ssp; $j++) {
my $id = $$srcIdCol[$j];
my $sc = $j+1;
info("**** Source ID $id ($sc of $num_ssp): Step 3 background");
if ( $X1_ima[$j] > 0.0 && $Y1_ima[$j] > 0.0 ) {
my $m = ($smallWindowMode && $sCCD[$j] == 1) ? 45 : 0;
while ($m < $MAX_M) {
$m = $m + 1;
info("[3] Source ID: $id, iteration $m: Running select_bkg(".$X1_ima[$j].",".$Y1_ima[$j].",".$R1_ima[$j].")");
my $result = &select_bkg($X1_ima[$j], $Y1_ima[$j], $R1_ima[$j], $m, $extrDetMask, $epicImage, $id);
$illum = $$result->{illum};
my $max_index = $$result->{index};
if ($illum > $good_area_limit) {
info("Source $id: CONVERGED @ radial step: $m azimuthal step: $max_index illumination: $illum");
$XB1_ima[$j] = $$result->{xback};
$YB1_ima[$j] = $$result->{yback};
$RB1_ima[$j] = $$result->{radback};
info("Source: $id: background solution XB1_ima, YB1_ima, RB1_ima: ".$XB1_ima[$j] . "," . $YB1_ima[$j] . "," . $RB1_ima[$j]);
last;
} elsif ($illum > $max_illum) {
$XB1_ima[$j] = $$result->{xback};
$YB1_ima[$j] = $$result->{yback};
$RB1_ima[$j] = $$result->{radback};
$max_illum = $illum;
}
}
}
$XB1_raw[$j] = ($XB1_ima[$j] - $crpix1) * 80.0 + $refx1;
$YB1_raw[$j] = ($YB1_ima[$j] - $crpix2) * 80.0 + $refy1;
$RB1_raw[$j] = $RB1_ima[$j] * 80.0;
$max_illum = 0;
}
info("==== Step 3 iteration complete");
info("==== Step 4: Start final extraction region generation");
my @regions;
for ($j = 0; $j < $num_ssp; $j++) {
my $id = $$srcIdCol[$j];
my $snum = $j+1;
info("* * * Writing region for source ID: $id ($snum of $num_ssp)");
my $sx = $X1_raw[$j];
my $sy = $Y1_raw[$j];
my $sr = $R1_raw[$j];
my $srcFITSRegionFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'Source extraction region file'
, format => 'FITS'
);
writeRegionRow($srcFITSRegionFile, "CIRCLE", $sx, $sy, $sr);
my $srcASCRegionFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'src ds9 region'
, format => 'ASCII'
);
my $rtn = writeASCIIRegion( $srcASCRegionFile
, $sx
, $sy
, $sr
, "green"
, "$id"
, 0 );
my $bx = $XB1_raw[$j];
my $by = $YB1_raw[$j];
my $br = $RB1_raw[$j];
my $bkgFITSRegionFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'Background extraction region file'
, format => 'FITS'
);
writeRegionRow($bkgFITSRegionFile, "CIRCLE", $bx, $by, $br);
my $bkgASCRegionFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'bkg ds9 region'
, format => 'ASCII'
);
my $rtn = writeASCIIRegion( $bkgASCRegionFile
, $bx
, $by
, $br
, "white"
, "$id"
, 0 );
my $k;
for ($k=0; $k < $num_sources; $k++) {
if ( $id == $$mlIdCol[$k] ) { next; }
my $cx = $allX_raw[$k];
my $cy = $allY_raw[$k];
my $cr = $allR_raw[$k];
my $sep = sqrt(($bx - $cx)**2 +
($by - $cy)**2);
if ($sep < ($br + $cr)) {
info("Contaminating source at $cx, $cy");
writeRegionRow($bkgFITSRegionFile, "!CIRCLE", $cx, $cy, $cr);
my $rtn = writeASCIIRegion( $bkgASCRegionFile
, $cx
, $cy
, $cr
, "red"
, "$k"
, 1 );
}
}
push (@regions, { 'id' => $id,
'source' => $srcFITSRegionFile,
'background' => $bkgFITSRegionFile
});
}
info("==== Step 4: Final extraction region generation complete");
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);
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");
my $id=$$srcIdCol[$i];
my $src_ra=$$srcRaCol[$i];
my $src_dec=$$srcDecCol[$i];
my $row_in_filteredML=$mlId_to_rowNum{ $mlSrcIdColInSSP->[$i] };
my $srcRegionSet = findFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Source extraction region file'
, format => 'FITS'
);
my $bkgRegionSet = findFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Background extraction region file'
, format => 'FITS'
);
my $echeckregion_output = newFile(class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'echeckregion output'
, format => 'ASCII'
);
my $echeckregion_status = ModuleUtil::isSourceOnChip(
"region($srcRegionSet,X,Y)"
, $forTimeFilteredEvList
, $echeckregion_output
);
if (!$echeckregion_status)
{
info("Source lies outside detector window, skipping this source for this exposure.");
next;
} else {
info("Source seems to be within detector window. Continue making products.");
}
my $srcRateSet = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Uncorrected source time series'
);
my $bgdRateSet = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Uncorrected background time series'
);
doCommand(
'evselect'
, table => "${forTimeFilteredEvList}:EVENTS"
, writedss => 'yes'
, updateexposure => 'yes'
, expression => "region($srcRegionSet,X,Y)"
, energycolumn => 'PI'
, withrateset => 'yes'
, rateset => $srcRateSet
, timebinsize => $$timeBinSizeCol[$i]
, maketimecolumn => 'yes'
, makeratecolumn => 'yes'
)
or return exception();
doCommand(
'evselect'
, table => "${forTimeFilteredEvList}:EVENTS"
, writedss => 'yes'
, updateexposure => 'yes'
, expression => "region($bkgRegionSet,X,Y)"
, energycolumn => 'PI'
, withrateset => 'yes'
, rateset => $bgdRateSet
, timebinsize => $$timeBinSizeCol[$i]
, maketimecolumn => 'yes'
, makeratecolumn => 'yes'
)
or return exception();
my @rate = readFITSColumn(file => $srcRateSet
, extension => 'RATE'
, column => 'RATE'
);
unless(grep { $_ } @rate)
{
info("Lightcurve contains no useful bins, skipping this source for this exposure");
next;
}
my $bkgSubtractedRateSet = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Exposure-corrected time series'
);
my $tempTS = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Temporary time series'
);
my $tempEventList = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Temporary event list'
);
my $lccorr_succeeded = doCommand(
'lccorr_pcms'
, srctsset => $srcRateSet
, eventset => $forTimeFilteredEvList
, withbkgtsset => 'yes'
, bkgtsset => $bgdRateSet
, subtractbkg => 'yes'
, treatvignet => 'no'
, treatfiltertrans => 'no'
, treatquantumeff => 'no'
, treatcosmicrays => 'yes'
, treatdeadtime => 'yes'
, treatgti => 'yes'
, treatalias => 'no'
, applycorrections => 'yes'
, outset => $bkgSubtractedRateSet
, srcweightstyle => 'deduce'
, srcra => $src_ra
, srcdec => $src_dec
, srcdeducestyle => 'events'
, patternlo => 0
, patternhi => $patternHi
, flagmask => $eventFlagMask
, pilos => 200
, pihis => 12000
, withgtis => 'yes'
, gtitabsnode0 => [@gtiTabList]
, rawylo => $rawylo
, srcregiontab => "$srcRegionSet:REGION"
, bkgregiontab => "$bkgRegionSet:REGION"
, tempset => $tempTS
, tempeventset => $tempEventList
);
foreach my $rebin (qw(2 4)) {
info("Generating x${rebin} rebinned source timeseries...");
my $srcRateSetN = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => "Uncorrected source time series (x${rebin} rebin)"
);
my $bgdRateSetN = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => "Uncorrected background time series (x${rebin} rebin)"
);
doCommand(
'evselect'
, table => "${forTimeFilteredEvList}:EVENTS"
, writedss => 'yes'
, updateexposure => 'yes'
, expression => "region($srcRegionSet,X,Y)"
, energycolumn => 'PI'
, withrateset => 'yes'
, rateset => $srcRateSetN
, timebinsize => $rebin * $$timeBinSizeCol[$i]
, maketimecolumn => 'yes'
, makeratecolumn => 'yes'
)
or return exception();
doCommand(
'evselect'
, table => "${forTimeFilteredEvList}:EVENTS"
, writedss => 'yes'
, updateexposure => 'yes'
, expression => "region($bkgRegionSet,X,Y)"
, energycolumn => 'PI'
, withrateset => 'yes'
, rateset => $bgdRateSetN
, timebinsize => $rebin * $$timeBinSizeCol[$i]
, maketimecolumn => 'yes'
, makeratecolumn => 'yes'
)
or return exception();
my $num_bins = numberFITSRows(file => $srcRateSetN
, extension => 'RATE');
unless($num_bins >= 5)
{
info("x${rebin} rebinned lightcurve contains insufficient bins ($num_bins), skipping...");
next;
}
info("Generating x${rebin} rebinned lightcurves...");
my $bkgSubtractedRateSetN = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => "Exposure-corrected time series (x${rebin} rebin)"
);
my $tempTSN = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => "Temporary time series (x${rebin} rebin)"
);
my $tempEventListN = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => "Temporary event list (x${rebin} rebin)"
);
my $lccorr_succeededN = doCommand(
'lccorr_pcms'
, srctsset => $srcRateSetN
, eventset => $forTimeFilteredEvList
, withbkgtsset => 'yes'
, bkgtsset => $bgdRateSetN
, subtractbkg => 'yes'
, treatvignet => 'no'
, treatfiltertrans => 'no'
, treatquantumeff => 'no'
, treatcosmicrays => 'yes'
, treatdeadtime => 'yes'
, treatgti => 'yes'
, treatalias => 'no'
, applycorrections => 'yes'
, outset => $bkgSubtractedRateSetN
, srcweightstyle => 'deduce'
, srcra => $src_ra
, srcdec => $src_dec
, srcdeducestyle => 'events'
, patternlo => 0
, patternhi => $patternHi
, flagmask => $eventFlagMask
, pilos => 200
, pihis => 12000
, withgtis => 'yes'
, gtitabsnode0 => [@gtiTabList]
, rawylo => $rawylo
, srcregiontab => "$srcRegionSet:REGION"
, bkgregiontab => "$bkgRegionSet:REGION"
, tempset => $tempTSN
, tempeventset => $tempEventListN
);
};
my $process_ts_further = 0;
if (!$lccorr_succeeded || !fileExists(file => $bkgSubtractedRateSet))
{
info("Lccorr failed for source $id");
}
else
{
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 EP_FLAG PN_FLAG M1_FLAG M2_FLAG EXTRRAD)]
,attributetype => [qw(integer real real real real string string string string real)]
,attributecomment => [
'EPIC source list number'
,'RA of source'
,'DEC of source'
,'Corrected RA of source'
,'Corrected DEC of source'
,'EP automatic source flag'
,'PN automatic source flag'
,'M1 automatic source flag'
,'M2 automatic source flag'
,'Source extraction radius (arcsec)'
]
,integervalue => $id
,realvalue => [ $$srcRaCol[$i] , $$srcDecCol[$i] , $$srcCorrRaCol[$i] , $$srcCorrDecCol[$i] , $source_radius[$id] ]
,stringvalue => [ $$epFlagCol[$i] , $$pnFlagCol[$i], $$m1FlagCol[$i], $$m2FlagCol[$i] ]
, set => "$bkgSubtractedRateSet:RATE"
) or exception();
my $mergedGtis = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Merged GTIs'
);
my @gtiTableList = ("${bkgSubtractedRateSet}:SRC_GTIS"
, "${bkgSubtractedRateSet}:BKG_GTIS");
if ($flareGtiExists) {push @gtiTableList, "${flareGtifile}:STDGTI";}
doCommand(
'gtimerge'
, tables => [@gtiTableList]
, withgtitable => 'yes'
, gtitable => $mergedGtis
, mergemode => 'and'
)
or return exception();
$alignedGtis = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Aligned merged GTIs'
);
doCommand(
'gtialign'
, style => 'generic'
, ingtitable => "$mergedGtis:STDGTI"
, outgtitable => "$alignedGtis:STDGTI"
, tstable => "${bkgSubtractedRateSet}:RATE"
)
or return exception();
my $tempCorrectedRateSet = newFile(
class => 'intermediate'
, instrument => thisInstrument
, band => 8
, src_num => $id
, exp_id => $exp_id
, content => 'Source timeseries'
);
doCommand( 'addattribute'
, attributename => [qw(AVRATE CHISQUAR CHI2PROB N_POINTS)]
, attributetype => [qw(real real real integer)]
, attributecomment => [
'Average rate in light curve'
, 'Chi squared value'
, 'Probability that CHISQUAR arises by chance'
, 'The number of data points'
]
, realvalue => [ -1.0, -1.0, -1.0 ]
, integervalue => [ 0 ]
, set => "$bkgSubtractedRateSet:RATE"
) or exception();
my $ekstest_output = doCommand(
'ekstest'
, set => $bkgSubtractedRateSet
, gtis => 'yes'
, gtiset => $alignedGtis
, chisquaretest => 'yes'
, fracvartest => 'no'
, screen => 'yes'
, newoutset => 'yes'
, outset => $tempCorrectedRateSet
)
or return exception();
my $ekstest_output = doCommand(
'ekstest'
, set => $tempCorrectedRateSet
, gtis => 'yes'
, gtiset => $alignedGtis
, chisquaretest => 'no'
, fracvartest => 'yes'
, screen => 'yes'
, newoutset => 'yes'
, outset => $tempCorrectedRateSet
)
or return exception();
my $fftPsFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'Time series FFT plot'
, format => 'PS'
);
my $efftplot_succeeded = doCommand(
'efftplot'
, infile => $tempCorrectedRateSet
, xscale => 'log'
, yscale => 'lin'
, gtis => 'yes'
, gtiset => $alignedGtis
, plotdevice => '/vcps'
, outfile => $fftPsFile
);
if (!$efftplot_succeeded || !fileExists(file => $fftPsFile))
{
info("efftplot failed for source $id");
}
}
my $srcSelExpression = "(region($srcRegionSet,X,Y))";
my $bkgSelExpression = "(region($bkgRegionSet,X,Y))";
if ($flareGtiExists) {
$srcSelExpression .= " && (GTI(${flareGtifile}:STDGTI,TIME))";
$bkgSelExpression .= " && (GTI(${flareGtifile}:STDGTI,TIME))";
}
my $tempSrcSpectrum = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Source spectrum'
);
my $tempBkgSpectrum = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'Source background spectrum'
);
my $srcSpectrum = newFile(
class => 'product'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'EPIC source spectrum'
) or exception();
my $bkgSpectrum = newFile(
class => 'product'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'EPIC source background spectrum'
);
my $arf = newFile(
class => 'intermediate'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'EPIC ancillary response function'
);
my $arfProd = newFile(
class => 'product'
, instrument => thisInstrument
, src_num => $id
, exp_id => $exp_id
, content => 'EPIC ancillary response function'
);
doCommand(
'especget'
, table => $forSpecFilteredEvList
, srcexp => $srcSelExpression
, backexp => $bkgSelExpression
, withfilestem => 'no'
, srcspecset => $tempSrcSpectrum
, bckspecset => $tempBkgSpectrum
, srcarfset => $arf
, witharfset => 'yes'
, withrmfset => 'yes'
, useodfatt => 'no'
, withbadpixcorr => 'yes'
, extendedsource => 'no'
)
or return exception();
doCommand( 'addattribute'
,attributename => [qw(SRCNUM SRC_RA SRC_DEC SRC_CRA SRC_CDEC EP_FLAG PN_FLAG M1_FLAG M2_FLAG EXTRRAD)]
,attributetype => [qw(integer real real real real string string string string real)]
,attributecomment => [
'EPIC source list number'
,'RA of source'
,'DEC of source'
,'Corrected RA of source'
,'Corrected DEC of source'
,'EP automatic source flag'
,'PN automatic source flag'
,'M1 automatic source flag'
,'M2 automatic source flag'
,'Source extraction radius (arcsec)'
]
,integervalue => $id
,realvalue => [ $$srcRaCol[$i] , $$srcDecCol[$i] , $$srcCorrRaCol[$i] , $$srcCorrDecCol[$i] , $source_radius[$id] ]
,stringvalue => [ $$epFlagCol[$i] , $$pnFlagCol[$i] , $$m1FlagCol[$i], $$m2FlagCol[$i] ]
, set => ["$tempSrcSpectrum:SPECTRUM" , "$tempBkgSpectrum:SPECTRUM", "$arf:SPECRESP" ]
) or exception();
if ( !readFITSKeyword( file => $tempSrcSpectrum, extension => 2, keyword => "EXPOSURE" ) ) {
info("No EXPOSURE keyword found in $tempSrcSpectrum. Skipping further processing for this source, " .
"and copying to product spectrum $srcSpectrum");
copyFile(source => "$tempSrcSpectrum", destination => $srcSpectrum);
next;
}
my $cannedRmfFileName = readFITSKeyword(
file => $tempSrcSpectrum
, extension => 'SPECTRUM'
, keyword => 'RESPFILE'
);
info("Canned RMF is $cannedRmfFileName");
my $cannedRmfPathName = findRmfFile( file => $cannedRmfFileName );
unless ( $cannedRmfPathName )
{
exception('MISSINGRMF','Unable to locate requested rmf : ', $cannedRmfFileName );
}
info("Absolute path to canned RMF is $cannedRmfPathName");
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";
}
doCommand(
'grppha'
, infile => $tempSrcSpectrum
, outfile => "\!$tempSrcSpectrum"
, chatter => 0
, comm => $grpphaCommExpr
)
or return exception();
my $specQuality = readFITSColumn( file => $tempSrcSpectrum
, extension => 'SPECTRUM'
, column => 'QUALITY'
);
my $specCounts = readFITSColumn( file => $tempSrcSpectrum
, extension => 'SPECTRUM'
, column => 'COUNTS'
);
my $specNumRows = numberFITSRows( file => $tempSrcSpectrum
, extension => 'SPECTRUM'
);
my $specTotalCounts = 0;
for ($i = 0; $i < $specNumRows; $i++) {
if ( $$specQuality[$i] == 0 ) {
$specTotalCounts += $$specCounts[$i];
}
}
info("Source $id: Got $specTotalCounts total good counts at QUALITY==0");
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"
);
my $spectrumPdfFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'EPIC source spectrum plot'
, format => 'PDF'
);
my $xpecCommandFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'XSPEC command file'
, format => 'ASCII'
);
if ( $specTotalCounts < 100 ) {
info("Source $id: Fewer than 100 good counts in source spectrum, skipping product generation");
} else {
info("Source $id: More than 100 good counts in source spectrum, generating products");
copyFile(
source => $arf
, destination => $arfProd
);
copyFile(
source => $tempSrcSpectrum
, destination => $srcSpectrum
);
copyFile(
source => $tempBkgSpectrum
, destination => $bkgSpectrum
);
my $tempCorrectedRateSet = findFile(
class => 'intermediate'
, instrument => thisInstrument
, band => 8
, src_num => $id
, exp_id => $exp_id
, content => 'Source timeseries'
);
my $correctedRateSet = newFile(
class => 'product'
, instrument => thisInstrument
, band => 8
, src_num => $id
, exp_id => $exp_id
, content => 'EPIC source timeseries'
);
if (fileExists(file => $tempCorrectedRateSet)) {
copyFile(
source => $tempCorrectedRateSet
, destination => $correctedRateSet
);
if (fileExists( file => $srcRegionSet ) ) {
doCommand(
'ftappend'
, infile => "${srcRegionSet}[REGION]"
, outfile => $correctedRateSet
) or return exception();
doCommand(
'fparkey'
, value => "REGION_SRC"
, fitsfile => "${correctedRateSet}[REGION]"
, keyword => "EXTNAME"
) or return exception();
};
if (fileExists( file => $bkgRegionSet) ) {
doCommand(
'ftappend'
, infile => "${bkgRegionSet}[REGION]"
, outfile => $correctedRateSet
) or return exception();
doCommand(
'fparkey'
, value => "REGION_BKG"
, fitsfile => "${correctedRateSet}[REGION]"
, keyword => "EXTNAME"
) or return exception();
};
my $tsPsFile = newFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'Source timeseries plot'
, format => 'PS'
);
doCommand(
'elcplot'
, set => $correctedRateSet
, usegtiset => 'yes'
, useflareset => 'no'
, gtiset => $alignedGtis
, offset => 'yes'
, offsetstyle => 'auto'
, forcestart => 'yes'
, binsize => 1
, plotdevice => '/vcps'
, plotfile => $tsPsFile
)
or return exception();
my $tsPdfFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, band => 8
, src_num => $id
, content => 'EPIC source timeseries plot'
, format => 'PDF'
);
PStoPDF(
source => $tsPsFile
, destination => $tsPdfFile
) or return exception();
my $srcASCRegionFile = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'src ds9 region'
, format => 'ASCII'
);
my $bkgASCRegionFile = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'bkg ds9 region'
, format => 'ASCII'
);
my $ds9RegionFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'EPIC source ds9 region'
, format => 'ASCII'
);
catASCIIFiles(
infiles => [$srcASCRegionFile, $bkgASCRegionFile, $allSrcRegionFile]
, outfile => $ds9RegionFile
) or return exception();
};
my $fftPsFile = findFile(
class => 'intermediate'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'Time series FFT plot'
, format => 'PS'
);
my $fftPdfFile = newFile(
class => 'product'
, instrument => thisInstrument
, exp_id => $exp_id
, src_num => $id
, content => 'EPIC source FFT plot'
, format => 'PDF'
);
if(fileExists(file => $fftPsFile)) {
PStoPDF(
source => $fftPsFile
, destination => $fftPdfFile
) or return exception();
}
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 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()
;
my $ancrfile = readFITSKeyword(
file => $srcSpectrum
, extension => 'SPECTRUM'
, keyword => 'ANCRFILE'
);
if ( $ancrfile )
{
my $ancrFileName = $arfProd;
$ancrFileName =~ s/\.FIT/\.FTZ/g;
doCommand('addattribute'
,set => "$srcSpectrum:SPECTRUM"
,attributename => 'ANCRFILE'
,attributetype => 'string'
,stringvalue => $ancrFileName
,attributecomment => "\"\\\"ancillary response\\\"\""
) or exception();
}
my $ancrFileName = $arfProd;
$ancrFileName =~ s/\.FIT/\.FTZ/g;
$ancrFileName =~ s/product\///g;
$bkgSpectrumFileName =~ s/\.FIT/\.FTZ/g;
$grpphaCommExpr = "chkey ANCRFILE $ancrFileName"."&chkey RESPFILE "
."$cannedRmfFileName"."&chkey BACKFILE $bkgSpectrumFileName"."&exit";
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.");
}
}
}
return success();
}
sub coordConv
{
info("--- Coordinate conversion subroutine ---");
my ($image, $ra, $decl) = @_;
if ($image eq '') { Exception->error('Module', "coordConv: Image not defined") };
if ($ra !~ /^\d+(\.\d*)?/) { Exception->error('Module', "coordConv: RA $ra not valid") };
if ($decl !~ /^-?\d+(\.\d*)?/) { Exception->error('Module', "coordConv: Dec $decl not valid") };
my $cmd = "ecoordconv imageset=$image x=$ra y=$decl coordtype=eqpos withcoords=yes";
info("CMD: $cmd");
my $output = `$cmd`;
info("$output");
my ($x,$y);
($x, $y) = ($output =~ /X:\s+Y:\s+(\d+)\s+(\d+)/);
my $sccd;
( $sccd ) = ($output =~ /centred on CCD: (\d+)/);
return ($x, $y, $sccd);
}
sub writeRegionRow
{
my ($file, $shape, $x, $y, $r) = @_;
my ($data, $ff, $nr, $status);
my $exp_id = exposureID( instrument => thisInstrument
, stream => thisStream
);
my @shape = ( $shape );
my @x = ( $x );
my @y = ( $y );
my @r = ( $r );
if (!fileExists( file => $file )) {
info("writeRegionRow(): Region file does not exist: creating it");
my $extrRegion = findFile(
class => "intermediate"
, instrument => thisInstrument
, exp_id => $exp_id
, content => 'extraction region file'
, format => 'FITS'
) or return exception();
doCommand(
'fselect'
, infile => $extrRegion
, outfile => $file
, expr => "SHAPE .eq. 'FOO'"
);
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'HDUVERS'
, attributetype => 'string'
, stringvalue => '1.0.0'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'HDUCLASS'
, attributetype => 'string'
, stringvalue => 'ASC'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'HDUCLAS1'
, attributetype => 'string'
, stringvalue => 'REGION'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'HDUCLAS2'
, attributetype => 'string'
, stringvalue => 'STANDARD'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'MTYPE1'
, attributetype => 'string'
, stringvalue => 'pos'
) or exception();
doCommand(
'addattribute'
, set => "$file:REGION"
, attributename => 'MFORM1'
, attributetype => 'string'
, stringvalue => 'X,Y'
) or exception();
}
fits_open_file($ff, $file, READWRITE, $status);
FITSException->error($status,'Error reading ', file => $file) if ($status);
$ff->movnam_hdu(ANY_HDU,'REGION', 0, $status);
FITSException->error($status,'Error accessing extension ', file => $file, extension => 'REGION') if ($status);
$ff->get_num_rows($nr, $status);
FITSException->error($status,'Error reading num_rows ', file => $file, extension => 'REGION') if ($status);
info("DEBUG: REGION file $file: writing $shape, $x, $y, $r firstrow $nr");
$ff->write_col(TSTRING, 1, $nr+1, 1, 1, \@shape, $status);
FITSException->error($status,'Error writing SHAPE ', file => $file, extension => 'REGION') if ($status);
$ff->write_col(TFLOAT, 2, $nr+1, 1, 1, \@x, $status);
FITSException->error($status,'Error writing X ', file => $file, extension => 'REGION') if ($status);
$ff->write_col(TFLOAT, 3, $nr+1, 1, 1, \@y, $status);
FITSException->error($status,'Error writing Y ', file => $file, extension => 'REGION') if ($status);
$ff->write_col(TFLOAT, 4, $nr+1, 1, 1, \@r, $status);
FITSException->error($status,'Error writing R ', file => $file, extension => 'REGION') if ($status);
$ff->close_file($status);
}
sub writeASCIIRegion
{
my ( $name, $x, $y, $r, $colour, $label, $append ) = @_;
if ($name eq '') { Exception -> error( 'Module', "writeAsciiRegion: filename not defined") };
if ($x !~ /^\d+(\.\d*)?/) { Exception -> error( 'Module', "writeAsciiRegion: detector X $x not valid") };
if ($y !~ /^\d+(\.\d*)?/) { Exception -> error( 'Module', "writeAsciiRegion: detector Y $y not valid") };
if ($r !~ /^\d+(\.\d*)?/) { Exception -> error( 'Module', "writeAsciiRegion: region radius R $r not valid") };
my $labstring = "";
if ($label ne "") { $labstring = " text={${label}}"; }
if ($colour eq '') { $colour = 'white'; }
my $colstring = "# color=".$colour;
my @t = ();
if ($append == 0) {
push (@t, "# Region file format: DS9 version 3.0\n");
push (@t, "global font = \"helvetica 10 normal\" edit=0 move=0 delete=1 include=1 fixed=0\n");
push (@t, "detector;circle($x, $y, $r) ${colstring}${labstring}\n");
} else {
push (@t, "detector;-circle($x, $y, $r) ${colstring}${labstring}\n");
}
writeASCIIFile( name => $name
, text => \@t
, append => $append
);
my $rtn = 1;
return $rtn;
}
1;