package Subs::BATSurvey;
##############################################################################
#
# DESCRIPTION: This subroutine does calibration of BAT event data.
# DESCRIPTION: It runs bateconvert to do energy scale calibrations,
# DESCRIPTION: and then runs batmaskwtevt to calculate mask weights for
# DESCRIPTION: each event which correspond to the burst position.
#
# HISTORY:
# HISTORY: $Log: BATSurvey.pm,v $
# HISTORY: Revision 1.12 2013/08/19 20:29:24 apsop
# HISTORY: Only log fetch_from_repository result if fetch succeeded, to
# HISTORY: avoid a Perl message about using an undefined variable. (If
# HISTORY: unsuccessful, fetch_from_repository already logs an error.)
# HISTORY:
# HISTORY: Revision 1.11 2013/05/30 07:54:54 apsop
# HISTORY: Log fetches from BAT repository
# HISTORY:
# HISTORY: Revision 1.10 2004/08/16 15:23:09 apsop
# HISTORY: Turn on writing of HISTORY keywords.
# HISTORY:
# HISTORY: Revision 1.9 2004/05/06 20:02:34 dah
# HISTORY: Add version number back into the header comments.
# HISTORY:
# HISTORY: Revision 1.8 2004/04/28 13:48:37 dah
# HISTORY: Fix attitude file instrument type.
# HISTORY:
# HISTORY: Revision 1.7 2004/04/16 20:21:18 dah
# HISTORY: Begin using embedded history records
# HISTORY:
#
# VERSION: 0.0
#
#
##############################################################################
use Subs::Sub;
@ISA = ("Subs::Sub");
use strict;
sub new {
my $proto=shift;
my $self=$proto->SUPER::new();
$self->{DESCRIPTION}="BAT Survey Processing";
return $self;
}
##################
# METHODS:
##################
sub body {
my $self=shift;
my $log =$self->log();
my $filename=$self->filename();
my $procpar =$self->procpar();
my $jobpar =$self->jobpar();
##################################################
# calibrate the raw DPH files
#################################################
foreach my $raw ($filename->get("rawdph", "b", "sv*", "*") ) {
#####################################################
# the calibrated files are given the same name with
# a different mode. These are temporary files which
# we will remove later.
#####################################################
my ($inst, $mode, $index) = $filename->parse($raw, "rawdph");
my $calib = $filename->get("rawdph", $inst, "calb".$mode, $index);
$log->entry("Skipping calibration of $raw since we have no tool for it");
system("cp $raw $calib");
}
#####################################################
# now merge the calibrated DPHs
# this makes a bunch of .dph files
#####################################################
$self->merge_dph();
################################################
# now I can delete the calibrated unmerged DPHs
################################################
foreach my $calib ($filename->get("rawdph", "b", "calb*", "*") ) {
$log->entry("Deleting $calib");
unlink $calib;
}
##############################################################
# Make sure we have some DPHs - otherwise, what's the point?
##############################################################
my @dphs = $filename->get("dph", "b", "sv*", "*");
unless(@dphs) {
$log->entry("No DPH survey data");
return;
}
######################################################
# create a summed DPH for each DPH
# first we set up the summing FTOOL
# for now we use fsumrows. This will be replaced
# by something better later
#####################################################
my $summer = Util::HEAdas->new("batsumdph")
->params({rows => "-",
clobber => "yes",
chatter => 5,
history => "yes"});
################################
# loop over each dph and sum it
################################
foreach my $dph ($filename->get("dph", "b", "sv*", "*")) {
my ($inst, $mode, $index) = $filename->parse($dph, "dph");
my $summed = $filename->get("dph", $inst, "totl", $index);
unlink $summed;
$log->entry("Summing $dph to make $summed");
###############################
# run the summing tool
###############################
$summer->params({infile=>$dph,
outfile=>$summed})
->run();
} # end of loop to sum DPHs
#####################################################
# now create a summed DPI from the summed DPH
# we create a temporary DPI for each summed DPH
# and then sum the DPIs. Usually there will only be
# one DPH
#####################################################
my $dph2dpi = Util::HEAdas->new("batdph2dpi")
->params({rows => 1,
levels => "-",
clobber => "yes",
chatter => 5,
history => "yes"});
foreach my $summed ($filename->get("dph", "b", "totl", "*")) {
#######################################################
# make a partial DPI file which will get deleted later
#######################################################
my ($inst, $mode, $index) = $filename->parse($summed, "dph");
my $dpi = $filename->get("dpi", $inst, "part", $index);
unlink $dpi;
$log->entry("Converting $summed to $dpi");
$dph2dpi->params({infile=>$summed,
outfile=>$dpi})
->run();
} # end of loop extracting DPIs
######################################################
# now we merge the DPIS
######################################################
my @dpis = $filename->get("dpi", "b", "part", "*");
if(@dpis==1) {
##########################################
# only one DPI, so rename just rename it
# this is the normal case
##########################################
my $dpi = $dpis[0];
my ($inst, $mode, $index) = $filename->parse($dpi, "dpi");
my $summed = $filename->get("dpi", $inst, "totl", "001");
$log->entry("One DPI $dpi - renaming it to $summed");
rename $dpis[0], $summed;
} elsif(@dpis>1) {
#################################################
# multiple DPIs - we need to add them
# I'm tired and this isn't sposed happen, so I'm
# not going to implement it quite yet
##################################################
$log->error(1, "Multiple DPI merge not implemented");
}
#################################################
# now do all that we need to do to the DPI file
#################################################
$self->process_dpi();
} # end of body method
################################################################################
# merge the dphraw files to form a single dph file
################################################################################
sub merge_dph {
my $self=shift;
my $log =$self->log();
my $filename=$self->filename();
my $procpar =$self->procpar();
my $jobpar =$self->jobpar();
#####################################
# get a list of all the raw DPH files
#####################################
my @raws = $filename->get("rawdph", "b", "calb*", "*");
unless(@raws) {
$log->entry("No Raw DPH files to merge");
return;
}
#####################################################################
# although it will rarely happen we need to guard against the case
# where different files have different EBOUNDS extensions
# So we make a hash key out of the E_MIN column of all the ebounds
# exztensions and collect files according to that key
#####################################################################
my %groups=();
foreach my $raw (@raws) {
my @e_min = Util::FITSfile->new($raw, "EBOUNDS")
->cols("E_MIN")
->table;
my $key = join "|", @e_min;
$key =~ s/\s*//g;
my $group = $groups{$key};
unless($group) {
$group = [];
$groups{$key} = $group;
}
push @{$group}, ($raw);
} # end of loop collecting groups
####################################################
# now loop over all the groups we collected above
####################################################
my $index=1;
foreach my $group (values %groups) {
my $dph = $filename->get("dph", "b", "surv", $index++);
$log->entry("Merging @{$group} into $dph");
my $list = Util::FITSlist->new(@{$group});
$list->extension("BAT_DPH");
if($list->count() == 1) {
##############################
# just one file, so copy it
##############################
my $file =$list->file(0);
$log->entry("Copying $file to $dph");
open IN, "<".$list->file(0);
open OUT, ">$dph";
while(<IN>) { print OUT $_; }
close IN;
close OUT;
} else {
################################
# multiple files, so merge them
################################
$list->merge($dph);
Util::FITSfile->new($list->file(0), "EBOUNDS")
->append_to($dph);
Util::GTIlist->new(@{$group})
->extension("GTI")
->merge_and_append_to($dph);
}
} # end of loop over groups
} # end of merge_dph method
################################################################################
# process the single DPI file. So far we extract hot pixels, make a sky
# image and plot the image.
################################################################################
sub process_dpi {
my $self=shift;
my $log =$self->log();
my $filename=$self->filename();
my $procpar =$self->procpar();
my $jobpar =$self->jobpar();
####################################################
# get the DPI file name. At this point there should
# be only one.
####################################################
my $dpi = $filename->get("dpi", "b", "totl", "001");
unless( -f $dpi) {
$log->entry("No DPI file");
return;
}
###########################################################
# get the detector enable/disable map from the repository
###########################################################
my $time = Util::FITSfile->new($dpi, 0)->keyword("TSTART");
$log->entry("BATSurvey.pm requesting from repository " .
"enabcal: $time");
my $enabled = $filename->fetch_from_repository('enabcal', 'b', '*', $time);
if ($enabled) {
$log->entry("fetched from repository " .
"enabcal: \$enabled=$enabled");
}
##################################################
# now make the hot pixel "master" map
##################################################
my $quality = $filename->get('qualcal', 'b');
$log->entry("Making master quality map $quality");
unlink $quality;
Util::HEAdas->new("bathotpix")
->params({infile=>$dpi,
outfile=>$quality,
detmask=>"$enabled\[1; FLAG(1)\]",
keepfract=>0.98,
guardfract=>0.25,
guardval=>5.0,
applygaps=>"yes",
mergemask=>"yes",
zerothresh=> 20.0,
clobber=>"yes",
chatter=>5,
history=>"yes"})
->run();
###############################################
# now make an image from the DPI
###############################################
my $image = $filename->get("skyimage", "b", "none", "001");
$log->entry("Creating sky image $image from $dpi");
unlink $image;
Util::HEAdas->new("batfftimage")
->params({infile => $dpi,
outfile => $image,
attitude => $filename->get('attitude', 's'),
aperture => $filename->fetch_cal("aperture", "b"),
detmask => $quality,
maskoffx => 0.0,
maskoffy => 0.00,
maskoffz => 0.00,
bat_z => 0.00,
rebalance => "yes",
oversampx => 2,
oversampy => 2,
pcodemap => "NO",
pcodethresh => 0.01,
corrections => "autocollim",
teldef => $filename->fetch_cal("teldef", "b"),
clobber => "yes",
append => "no",
chatter => 5,
history => "yes"})
->run();
############################################
# now use ximage to plot the image
############################################
my $plot = $filename->corresponding("skyimage", "skyplot", $image);
$log->entry("Plotting $image in $plot");
Util::Xanadu->new("ximage")
->script("read/fits $image",
"smooth",
"cpd /vcps",
"cey 2000",
"disp",
"grid",
"scale",
"exit" )
->run();
rename "pgplot.ps", $plot;
} # end of process_dpi method