package Subs::UvotImages;
##############################################################################
#
# DESCRIPTION: Extract transform and plot images
#
# HISTORY:
# HISTORY: $Log: UvotImages.pm,v $
# HISTORY: Revision 1.47 2005/04/06 15:43:42 apsop
# HISTORY: Change to using CALDB for cal parameters.
# HISTORY:
# HISTORY: Revision 1.46 2005/03/16 13:31:37 apsop
# HISTORY: Add in FILTER keyword to image primary header.
# HISTORY:
# HISTORY: Revision 1.45 2005/03/15 18:58:39 apsop
# HISTORY: Put basic and DATE* keywords into image file primary headers. Fix bug (?) which was not using modmap output for processing.
# HISTORY:
# HISTORY: Revision 1.44 2005/02/08 18:21:32 apsop
# HISTORY: Add new types for uvot filter specific images and exp maps. Replace Filename::corresponding method with file name hack for now. Fix bug which was preventing the production of raw images.
# HISTORY:
# HISTORY: Revision 1.43 2005/01/21 04:26:43 apsop
# HISTORY: Change to use file names stored in cat file rather than contructing them.
# HISTORY:
# HISTORY: Revision 1.42 2005/01/12 17:29:23 apsop
# HISTORY: Change weightfile param in uvotdetect to expfile param.
# HISTORY:
# HISTORY: Revision 1.41 2005/01/08 04:11:00 apsop
# HISTORY: Change test for existence of event list to be more forgiving.
# HISTORY:
# HISTORY: Revision 1.40 2005/01/07 20:34:21 apsop
# HISTORY: No longer necessary to use the CENTER method of swiftxform or a small
# HISTORY: value of ncell for uvotmodmap.
# HISTORY:
# HISTORY: Revision 1.39 2004/12/03 13:38:11 apsop
# HISTORY: Turn back on calling uvotmodmap.
# HISTORY:
# HISTORY: Revision 1.38 2004/10/13 01:38:31 apsop
# HISTORY: Disable uvotmodmap for the time being, and fix bug in getting plot file names.
# HISTORY:
# HISTORY: Revision 1.37 2004/10/12 16:26:36 apsop
# HISTORY: Filename changes and a few other tweaks to make module run in complete pipeline.
# HISTORY:
# HISTORY: Revision 1.36 2004/09/15 22:26:24 apsop
# HISTORY: Remove the binning factor from uvot event list file names
# HISTORY:
# HISTORY: Revision 1.35 2004/09/01 14:42:21 apsop
# HISTORY: Disable level 2/3 processing for now.
# HISTORY:
# HISTORY: Revision 1.34 2004/08/31 20:51:02 apsop
# HISTORY: Implemented UVOT level 2 and 3 products.
# HISTORY:
# HISTORY: Revision 1.33 2004/08/18 17:27:58 apsop
# HISTORY: Yet another fix for the timing keywords.
# HISTORY:
# HISTORY: Revision 1.32 2004/07/19 16:03:24 apsop
# HISTORY: Add version number entry to comment field.
# HISTORY:
# HISTORY: Revision 1.31 2004/07/12 13:46:09 apsop
# HISTORY: Fix problems with timing keywords
# HISTORY:
# HISTORY: Revision 1.30 2004/06/15 16:01:04 apsop
# HISTORY: Better trapping of problem uvot exposures.
# HISTORY:
# HISTORY: Revision 1.29 2004/06/02 16:11:19 apsop
# HISTORY: Change name of catalogue extension.
# HISTORY:
# HISTORY: Revision 1.28 2004/05/04 16:33:01 dah
# HISTORY: Bug fixes for case where no initial image file exists.
# HISTORY:# HISTORY: Revision 1.27 2004/04/30 16:16:40 dah
# HISTORY: Use new FITSfile::keywords() method to speed up image extraction
# HISTORY:
# HISTORY: Revision 1.26 2004/04/16 20:21:18 dah
# HISTORY: Begin using embedded history records
# HISTORY:
# HISTORY: Revision 1.25 2004/04/16 19:53:18 dah
# HISTORY: Begin using embedded history records
# HISTORY:
#
# VERSION: 0.0
#
##############################################################################
use Subs::Sub;
use Subs::Images;
use Util::Xanadu;
use Subs::UvotNames;
@ISA = ("Subs::Images");
use strict;
sub new {
my $proto=shift;
my $self=$proto->SUPER::new();
$self->{DESCRIPTION}="Extracting, merging and plotting images for UVOT";
return $self;
}
##################
# METHODS:
##################
sub body {
my $self=shift;
#####################################
# extract images from the event data
#####################################
$self->extract_images();
$self->process_raw();
$self->image_raw_to_det();
$self->image_raw_to_sky('uvot');
$self->create_exposure_maps();
# create summed images/exposure maps
$self->sum_images();
# detect sources and determine magnitudes
$self->detect_sources();
$self->plot_images();
} # end of body method
#############################################################################
# Extract images from filtered event data
#############################################################################
sub extract_images {
my $self=shift;
my $filterCodes = $Subs::UvotNames::filterCodes;
my $filterNames = $Subs::UvotNames::filterParams;
my @exp_keywords = ('TSTART', 'TSTOP', 'TELAPSE', 'ONTIME', 'LIVETIME', 'EXPOSURE', 'DEADC');
my @ds_keywords = ('DSUNI1', 'DSVAL1', 'DSTYP1', 'DSREF2', 'DSUNI2', 'DSVAL2', 'DSTYP2');
my $log =$self->log();
my $filename=$self->filename();
my $procpar =$self->procpar();
my $jobpar =$self->jobpar();
$log->entry("Extracting images from uvot event data as required");
my $catFile = $filename->get('hk', 'uvot', 'ct', '*');
unless( -f $catFile ){
$log->entry("No catalogue file, so cannot extract images.");
return;
}
my $catFits = Util::FITSfile->new($catFile, 'EXPCATALOG');
unless ($catFits->keyword('EXTNAME')) {
$log->error(2, "File $catFile does not have an EXPCATALOG extension.");
return;
}
my @filterId = $catFits->cols('FILTER')->table();
my @gtiName = $catFits->cols('GTINAME')->table();
my @expRef = $catFits->cols('EXPREF')->table();
my @efile = $catFits->cols('EFILEREF')->table();
my @ifile = $catFits->cols('IFILEREF')->table();
my @binning = $catFits->cols('BINLVL')->table();
my @hasImage = $catFits->cols('Image')->table();
my @hasEvents = $catFits->cols('Event')->table();
my @startTime = $catFits->cols('ESTART')->table();
my @stopTime = $catFits->cols('ESTOP')->table();
my @x0 = $catFits->cols('EV_X0')->table();
my @y0 = $catFits->cols('EV_Y0')->table();
my @xsize = $catFits->cols('EV_XSIZ')->table();
my @ysize = $catFits->cols('EV_YSIZ')->table();
my @imx0 = $catFits->cols('IM_X0')->table();
my @imy0 = $catFits->cols('IM_Y0')->table();
my @imxsize = $catFits->cols('IM_XSIZ')->table();
my @imysize = $catFits->cols('IM_YSIZ')->table();
my @winx0 = $catFits->cols('WINDOWX0')->table();
my @winy0 = $catFits->cols('WINDOWY0')->table();
my @windx = $catFits->cols('WINDOWDX')->table();
my @windy = $catFits->cols('WINDOWDY')->table();
my @binx = $catFits->cols('BINX')->table();
my @biny = $catFits->cols('BINY')->table();
my $catRows = $catFits->nrows();
###########################
# fire up ftools
###########################
my $fextract = Util::Ftool->new('fextract');
my $tabedit = Util::HEAdas->new('ftedit');
my %needsSorting;
my $tmpFile = 'Uvot_raw_images.tmp';
for(my $exp=0; $exp<$catRows; $exp++){
if(!$binning[$exp] || $binning[$exp] eq 'INDEF'){
$log->entry("Exposure row $exp has an undefined or zero binning value." .
" Skipping extraction of this exposure.");
next;
}
next if $hasEvents[$exp] eq 'b0';
##################################################################
# Only make images for those exposures that don't already have one
# or if the event window is not completely contained in the image
# window
##################################################################
unless($hasImage[$exp] eq 'b0'){
my $extract = ($x0[$exp]-$xsize[$exp]/2) < ($imx0[$exp]-$imxsize[$exp]/2)
|| ($y0[$exp]-$ysize[$exp]/2) < ($imy0[$exp]-$imysize[$exp]/2)
|| ($x0[$exp]+$xsize[$exp]/2) > ($imx0[$exp]+$imxsize[$exp]/2)
|| ($y0[$exp]+$ysize[$exp]/2) > ($imy0[$exp]+$imysize[$exp]/2);
next unless $extract;
}
my $expref = $expRef[$exp];
$binning[$exp] =~ s/\s//g;
my $evtFile = $efile[$exp];
my $mode = $$filterCodes[$filterId[$exp]] . 'po';
my $imageFile = $filename->get('rawimage', 'uvot', $$filterCodes[$filterId[$exp]], 0);
my $extname = $$filterCodes[$filterId[$exp]] . sprintf('%08d', int($startTime[$exp])) . 'E';
my $evtFits = Util::FITSfile->new($evtFile, 'WINDOW');
my @windExpRef = $evtFits->cols('EXPREF')->table();
###################################
# Update catalogue and event files
###################################
$tabedit->params({infile => $catFile.'[1]',
row => $exp+1});
$tabedit->params({column => 'IFILEREF',
value => $imageFile})
->run();
$tabedit->params({column => 'IMENAME',
value => $extname})
->run();
$tabedit->params({column => 'Image',
value => 128})
->run();
for(my $irow=0; $irow<@windExpRef; $irow++){
if($windExpRef[$irow]==$expref){
$tabedit->params({infile => $evtFile.'[WINDOW]',
row => $irow+1,
column => 'IMAGEXT',
value => 1})
->run();
}
}
##############################
# Set up to extract the image
##############################
my $imageFits;
if( -f $imageFile ){
$needsSorting{$imageFile} = 1;
$imageFits = Util::FITSfile->new($imageFile, 0);
}else{
Util::Ftool->new('fimgcreate')
->params({bitpix => 8,
naxes => '0',
datafile => 'none',
outfile => $imageFile})
->run();
$evtFits->ext(0);
$imageFits = Util::FITSfile->new($imageFile, 0);
$imageFits->keyword('TELESCOP', $evtFits->keyword('TELESCOP'));
$imageFits->keyword('INSTRUME', $evtFits->keyword('INSTRUME'));
$imageFits->keyword('FILTER', $evtFits->keyword('FILTER'));
}
$evtFits->ext('EVENTS');
$evtFits->specs('[gtifilter(\'[' . $gtiName[$exp] . ']\')]' .
'[bin RAWX=' . int($winx0[$exp]) . ':' .
int($winx0[$exp]+$windx[$exp]-1) . ':' .
int($binx[$exp]) . ',' .
'RAWY=' . int($winy0[$exp]) . ':' .
int($winy0[$exp]+$windy[$exp]-1) . ':' .
int($biny[$exp]) . ']');
$log->entry("Extracting raw image from $evtFile, gti " . $gtiName[$exp] .' '. `date`);
$evtFits->append_to($imageFile);
$imageFits->ext($imageFits->nhdus()-1);
###############################################################
# Fix up the keywords. You'd think this would be easy ...
###############################################################
my %keywords = $imageFits->keywords();
my @keynames = keys %keywords;
$imageFits->begin_many_keywords();
$imageFits->keyword('EXTNAME', $extname);
my ($cdelt1, $cdelt2) = ($imageFits->keyword('CDELT1'), $imageFits->keyword('CDELT2'));
$imageFits->keyword('WINDOWX0', 0.5 + (0.5-$imageFits->keyword('CRPIX1'))*$cdelt1 +
$imageFits->keyword('CRVAL1') );
$imageFits->keyword('WINDOWY0', 0.5 + (0.5-$imageFits->keyword('CRPIX2'))*$cdelt1 +
$imageFits->keyword('CRVAL2') );
$imageFits->keyword('WINDOWDX', $imageFits->keyword('NAXIS1')*$cdelt1 );
$imageFits->keyword('WINDOWDY', $imageFits->keyword('NAXIS2')*$cdelt2 );
$imageFits->keyword('BINX', $cdelt1);
$imageFits->keyword('BINY', $cdelt2);
$imageFits->keyword('CMPUTHRS', -1000);
$imageFits->keyword('CMPOTHRS', 1000);
$imageFits->keyword('CMPCNTMN', 'F');
$imageFits->keyword('BLANK', -1);
$imageFits->keyword('TSTART', $startTime[$exp]);
$imageFits->keyword('TSTOP', $stopTime[$exp]);
$imageFits->keyword('TELAPSE', $stopTime[$exp]-$startTime[$exp] );
my $first = Util::Date->new($startTime[$exp]);
my $last = Util::Date->new($stopTime[$exp]);
$imageFits->keyword('DATE-OBS', $first->date().'T'.$first->time() );
$imageFits->keyword('DATE-END', $last->date().'T'.$last->time() );
foreach my $exp_key (@exp_keywords){
for(my $index=1; $index<=$catRows; $index++){
my $oldname = substr($exp_key, 0, 6) . $index;
$imageFits->keyword( "-$oldname", ' ' ) if grep /${oldname}/, @keynames;
if($index==$expref){
$imageFits->keyword( $exp_key, $evtFits->keyword($oldname) );
}
}
}
foreach my $ds_key (@ds_keywords){
for(my $index=1; $index<=$catRows; $index++){
my $oldname = $index==1 ? '' : $index;
$oldname .= $ds_key;
$imageFits->keyword( "-$oldname", ' ' ) if grep /${oldname}/, @keynames;
}
}
$imageFits->keyword( 'HDUCLAS1', 'IMAGE' );
$imageFits->keyword( 'HDUCLAS2', 'TOTAL' );
$imageFits->end_many_keywords();
my $livetime = $imageFits->keyword('ONTIME');
my $deadtime = $imageFits->keyword('DEADC');
$livetime = $livetime*$deadtime if $deadtime;
$imageFits->keyword( 'LIVETIME', $livetime);
$imageFits->keyword( 'EXPOSURE', $livetime);
}
#############################################################
# Sort the extensions to be in the right order, if necessary
#############################################################
my $fappend = Util::Ftool->new("fappend");
foreach my $sortFile (keys %needsSorting) {
$log->entry("Sorting image file $sortFile.");
unlink $tmpFile;
$fextract->params({'infile' => $sortFile . '[0]', 'outfile' => $tmpFile})
->run();
my @hduNames;
my $sortFits = Util::FITSfile->new($sortFile, 0);
my $nhdus = $sortFits->nhdus();
for(my $hdu=1; $hdu<$nhdus; $hdu++){
$sortFits->ext($hdu);
push @hduNames, $sortFits->keyword('EXTNAME');
}
my @sortedHduNames = sort @hduNames;
next if join('%',@sortedHduNames) eq join('%',@hduNames);
foreach my $extName (@sortedHduNames) {
next unless $extName;
$fappend->params({'infile' => $sortFile . "[$extName]",
'outfile' => $tmpFile})
->run();
}
rename $tmpFile, $sortFile;
}
unlink $tmpFile;
########################################################
# Calculate the timing keywords for the primary headers
########################################################
foreach my $rawFile ($filename->get('rawimage', 'uvot', '*', 0)){
my $rawFits = Util::FITSfile->new($rawFile);
my $nhdus = $rawFits->nhdus();
my ($tstart, $tstop) = (1.E10, 0);
for(my $hdu=1; $hdu<$nhdus; $hdu++){
$rawFits->ext($hdu);
my $start = $rawFits->keyword('TSTART');
my $stop = $rawFits->keyword('TSTOP');
$tstart = $start if $start < $tstart;
$tstop = $stop if $stop > $tstop;
}
my $first = Util::Date->new($tstart);
my $last = Util::Date->new($tstop);
$rawFits->ext(0);
# $rawFits->keyword('TSTART', $tstart);
# $rawFits->keyword('TSTOP', $tstop);
# $rawFits->keyword('TELAPSE', $tstop-$tstart);
$rawFits->keyword('DATE-OBS', $first->date().'T'.$first->time() );
$rawFits->keyword('DATE-END', $last->date().'T'.$last->time() );
}
} # end extract images method
###############################################################################
# uvot specific processing of the raw images
###############################################################################
sub process_raw {
my $self=shift;
my $log =$self->log();
my $filename=$self->filename();
my $procpar =$self->procpar();
my $jobpar =$self->jobpar();
$log->entry("Processing raw uvot images");
my $ubadpix = Util::HEAdas->new('uvotbadpix');
$ubadpix->params({'infile' => 'CALDB'});
my $umodmap = Util::HEAdas->new('uvotmodmap');
my $uflatfield = Util::HEAdas->new('uvotflatfield');
$uflatfield->params({flatfile => 'CALDB'});
$log->entry("setting {RA,DEC,PA}_PNT in all raw images");
my $ra = $jobpar->read("ra");
my $dec = $jobpar->read("dec");
my $roll = $jobpar->read("roll");
##############################
# iterate over the UVOT images
##############################
foreach my $rawFile ($filename->get('rawimage', 'uvot', '*', '*')) {
{
my $imageFits = Util::FITSfile->new($rawFile);
for my $i (1 .. $imageFits->nhdus - 1) {
$imageFits->ext($i);
$imageFits->keyword('RA_PNT', $ra);
$imageFits->keyword('DEC_PNT', $dec);
$imageFits->keyword('PA_PNT', $roll);
}
}
my $badFile = $filename->corresponding('rawimage', 'badimage', $rawFile);
my $corrFile = $filename->corresponding('rawimage', 'corrimage', $rawFile);
my $mod8File = 'mod.tmp';
$log->entry("running uvotbadpix on $rawFile");
$ubadpix->params({infile => $rawFile,
badpixlist => 'CALDB',
outfile => $badFile,
chatter => 3})
->run();
$log->entry("running uvotmodmap on $rawFile");
$umodmap->params({infile => $rawFile,
badpixfile => $badFile,
outfile => $mod8File,
nsig => 3,
ncell => 16, # TODO: keeps run times for 2048^2 images to minutes
# instead of hours
chatter => 3})
->run();
# uflatfield will go here
$uflatfield->params({
infile => $mod8File,
###### infile => $rawFile,
outfile => $corrFile})
->run();
unlink($mod8File);
}
}
###############################################################################
# convert raw coordinate images to detector coordinates
###############################################################################
sub image_raw_to_det {
my $self=shift;
my $log =$self->log();
my $filename=$self->filename();
my $procpar =$self->procpar();
my $jobpar =$self->jobpar();
$log->entry("Converting raw coordinate images to detector coordinates");
#########################################
# create the swiftxform tool
#########################################
my $swiftxform = Util::HEAdas->new('swiftxform');
$swiftxform->params({
teldeffile => 'CALDB',
to => 'DET',
attfile => 'NONE',
aberration => 'no',
seed => $procpar->read('seed'),
chatter => 4,
})
->is_script(1);
######################
# loop over _GRISM_ files
# since only they are converted to DET coordinates
######################
foreach my $rawFile ($filename->get('corrimage', 'uvot', 'g*', '*') ) {
#########################################
# determine the DET coordinate filename
#########################################
my $detFile = $filename->corresponding("corrimage", "detimage", $rawFile);
$log->entry("converting $rawFile to $detFile");
######################
# do the conversion
######################
$swiftxform->params({
infile => $rawFile,
outfile => $detFile,
})
->run();
} # end of loop over files
} # end of image_raw_to_det method
###############################################################################
# convert raw coordinate images to sky coordinates
###############################################################################
sub image_raw_to_sky {
my $self=shift;
my $log =$self->log();
my $filename=$self->filename();
my $procpar =$self->procpar();
my $jobpar =$self->jobpar();
$log->entry("Converting raw coordinate images to sky coordinates");
######################################################
# make sure there is an attitude file
######################################################
my $attitude = $filename->get('attitude', 's');
unless(-f $attitude) {
$log->entry("No attitude data available - can't make sky images");
return;
}
#########################################
# create the swiftxform tool
#########################################
my $swiftxform = Util::HEAdas->new('swiftxform')->is_script(1);
$swiftxform->params({
teldeffile => 'CALDB',
to => 'SKY',
attfile => $attitude,
ra => $jobpar->read('ra'),
dec => $jobpar->read('dec'),
method => 'AREA',
aberration => 'no',
seed => $procpar->read('seed'),
chatter => 4,
});
######################
# loop over files
######################
foreach my $rawFile ($filename->get('corrimage', 'uvot', '*', '*') ) {
#########################################
# determine the DET coordinate filename
#########################################
my $skyFile = $rawFile;
$skyFile =~ s/_cr\./_sk\./;
$log->entry("converting $rawFile to $skyFile");
######################
# do the conversion
######################
$swiftxform->params({
infile => $rawFile,
outfile => $skyFile,
})
->run();
} # end of loop over files
} # end of image_raw_to_sky method
###############################################################################
# generate exposure maps
###############################################################################
sub create_exposure_maps {
my $self=shift;
my $log = $self->log();
my $filename= $self->filename();
my $procpar = $self->procpar();
my $jobpar = $self->jobpar();
$log->entry("Creating exposure maps");
######################################################
# make sure there is an attitude file
######################################################
my $attitude = $filename->get('attitude', 's');
unless(-f $attitude) {
$log->entry("No attitude data available - can't make exposure maps");
return;
}
#########################################
# create the swiftxform tool
#########################################
my $uexpmap = Util::HEAdas->new('uvotexpmap');
$uexpmap->params({
attfile => $attitude,
teldeffile => 'CALDB',
method => 'MEANFOV', # TODO: is SHIFT/ADD in effect?
attdelta => 100, # TODO: make a parameter
aberration => 'no',
chatter => 4,
});
######################
# loop over files
######################
foreach my $skyFile ($filename->get('filterimg', 'uvot', '*', '*') ) {
# exclude grism files
next if $filename->is_grism($skyFile, 'skyimage');
#########################################
# determine the exposure map filename
#########################################
my $expFile = $skyFile;
$expFile =~ s/_sk\./_ex\./;
$log->entry("building exposure map $expFile for $skyFile");
my $badFile = $skyFile;
$badFile =~ s/_sk\./_bp\./;
######################
# do the conversion
######################
$uexpmap->params({
infile => $skyFile,
badpixfile => $badFile,
outfile => $expFile,
})
->is_script(1)
->run();
} # end of loop over files
} # end of create_exposure_maps method
sub sum_images {
my $self=shift;
my $log = $self->log();
my $filename= $self->filename();
my $procpar = $self->procpar();
my $jobpar = $self->jobpar();
$log->entry("Creating summed images");
################################
# Create tool for summing images
################################
my $usumexp = Util::HEAdas->new('uvotimsum');
my $tmpFile = 'sum.tmp';
########################################
# clean up any pre-existing summed files
########################################
foreach my $type (qw(skyimage expimage)) {
foreach my $imageFile ($filename->get($type, 'uvot', '', '*')) {
$log->entry("removing obsolete summed image $imageFile");
unlink($imageFile);
}
}
######################
# loop over files
######################
my %filter_to_sum = ('filterimg' => 'skyimage',
'filterexp' => 'expimage');
foreach my $type (keys %filter_to_sum) {
foreach my $imageFile ($filename->get($type, 'uvot', '*', '*')) {
my (undef, $filter_code) = $filename->parse($imageFile, $type);
$log->entry("Filter code for $imageFile is $filter_code.");
next if not $filter_code;
#########################################
# determine the summed image filename
# by removing the filter code from sw\d{9}uff(sk|ex).img
#########################################
my $sumFile = $filename->get($filter_to_sum{$type}, 'uvot', '', 0);
$log->entry("summing $imageFile to $sumFile");
my $filter = $filename->filter_name_for_code($filter_code);
######################
# do the summation
######################
$usumexp->params({
infile => $imageFile,
outfile => $tmpFile,
})
->is_script(1)
->run();
my $keys = "[col "
. " #EXTNAME='$filter';"
. " #FILTER='$filter';"
. ']';
if (not -e $sumFile) {
Util::HEAdas->new('ftcopy')
->params({
infile => $tmpFile . "[col #EXTNAME='$filter']",
outfile => $sumFile,
copyall => 'yes',
})
->run();
}
else {
# append this result to the sum file
Util::HEAdas->new('ftappend')
->params({
infile => "$tmpFile+1[col #EXTNAME='$filter']",
outfile => $sumFile,
})
->run();
}
unlink($tmpFile);
} # end of loop over files
} # end of loop over types
} # end of process_sky method
sub detect_sources {
my $self=shift;
my $log = $self->log();
my $filename= $self->filename();
my $procpar = $self->procpar();
my $jobpar = $self->jobpar();
$log->entry("Detecting sources");
########################################
# clean up any pre-existing source lists
########################################
foreach my $catFile ($filename->get('srclist', 'uvot', '', '*')) {
unlink($catFile);
}
###############################################################
# Create tools for detecting sources and determining magnitudes
###############################################################
my $udetect = Util::HEAdas->new('uvotdetect');
my $umag = Util::HEAdas->new('uvotmag');
$umag->params({
zerofile => 'CALDB',
coinfile => 'CALDB',
filter => 'default',
});
my $tmpFile = 'sources.tmp';
my $catFile = undef;
#############################
# loop over summed sky images
#############################
foreach my $skyFile ($filename->get('skyimage', 'uvot', '', '*')) {
$log->entry("detect: skyFile '$skyFile'");
my $expFile = $filename->corresponding('skyimage', 'expimage', $skyFile);
if (not $expFile or not -f $expFile) {
$log->entry("No exposure data for $skyFile - assuming constant exposure");
$expFile = 'NONE';
}
if (not $catFile) {
$catFile = $filename->corresponding('skyimage', 'srclist', $skyFile);
}
#########################################
# run uvotdetect on each extension
#########################################
my $fitsFile = Util::FITSfile->new($skyFile, 0);
my $nhdus = $fitsFile->nhdus();
for (my $hdu0 = 1; $hdu0 < $nhdus; ++$hdu0) {
$fitsFile->ext($hdu0);
my $filter = $fitsFile->keyword('EXTNAME');
my $weightFile = $expFile;
if ($expFile ne 'NONE') {
$weightFile .= "[$filter]";
}
if (-e $tmpFile) {
unlink($tmpFile);
}
$udetect->params({
infile => "$skyFile+$hdu0",
outfile => $tmpFile,
expfile => $weightFile,
threshold => 2, # TODO: make parameter
})
->is_script(1)
->run();
# run uvotmag
if (-e $tmpFile) {
$umag->params({
infile => $tmpFile,
})
->is_script(1)
->run();
}
if (not -e $catFile) {
Util::HEAdas->new('ftcopy')
->params({
infile => $tmpFile . "[col #EXTNAME='$filter']",
outfile => $catFile,
copyall => 'yes',
})
->run();
}
else {
# append this result to the sum file
Util::HEAdas->new('ftappend')
->params({
infile => "$tmpFile+1[col #EXTNAME='$filter']",
outfile => $catFile,
})
->run();
}
} # end of loop over extensions
} # end of loop over files
unlink($tmpFile) if -e $tmpFile;
} # end of detect_sources method
################################################################################
#
################################################################################
sub plot_images {
my $self = shift;
my $log =$self->log();
my $filename=$self->filename();
my $procpar =$self->procpar();
my $jobpar =$self->jobpar();
my $max_image_dimen=1024;
$log->entry("Plotting all images");
my $ximage = Util::Xanadu->new("ximage");
my $type;
foreach $type ("detimage", "skyimage") {
my $plot_type = $type;
$plot_type =~ s/image/plot/;
######################################
# loop over summed images of this type
######################################
my $file;
foreach my $file ($filename->get($type, 'uvot', '', '*')) {
my $fitsFile = Util::FITSfile->new($file, 0);
my $nhdus = $fitsFile->nhdus();
my $grid='';
if($type =~ /^sky/ ) { $grid = 'grid' }
##############################################
# for each extension (filter) create a graphic
##############################################
for(my $hdu=1; $hdu<$nhdus; $hdu++) {
my @commands;
push @commands, ('cpd /vgif', 'cey 2000');
$fitsFile->ext($hdu);
my $filter = $fitsFile->keyword('EXTNAME');
my $filter_code = $filename->filter_code_for_name($filter);
my $plotFile = $filename->get($plot_type, 'uvot', $filter_code, 0);
########################################
# determine if we need to bin the image
########################################
my $dimen = $fitsFile->keyword('NAXIS1');
my $bin=1;
for($bin=1; $max_image_dimen*$bin<$dimen; $bin++) {}
$log->entry("Plotting ${file}\[$hdu\] in $plotFile binned by $bin");
push @commands, ("read/fits/rebin=$bin \{${file}\[$hdu\]\}",
'smooth',
'disp',
$grid,
'scale');
push @commands, 'exit';
$ximage->script(@commands)->run();
rename 'pgplot.gif', $plotFile;
}
}
}
} # end of plot_images method
1;