package Subs::XrtGrbLc;
##############################################################################
#
# DESCRIPTION: Produce XRT GRB lightcurve, gif, and related files. Runs
# xrtgrblc for PC and WT mode events; xrtproducts for PD/LR and PD/PU mode
# events.
#
# HISTORY:
# HISTORY: $Log: XrtGrbLc.pm,v $
# HISTORY: Revision 1.16 2014/04/21 21:09:36 apsop
# HISTORY: Generate light curve, pha, and arf files for XRT PD/LR and
# HISTORY: PD/PU modes. New subs make_pd_mode_lc and concatLogs.
# HISTORY: Required some reorganization to allow processing PD modes
# HISTORY: potentially without PC or WT modes. write_standard_keyword
# HISTORY: now logs the file it's doing. Several minor and commentary
# HISTORY: improvements.
# HISTORY:
# HISTORY: Revision 1.15 2014/02/28 11:58:29 apsop
# HISTORY: Don't use the uat attitude file, at request of XRT Team.
# HISTORY: No longer sets ATTFLAG keywords (leaves that for SW0WrapUp).
# HISTORY: Check for final run using jobpar object rather than job_title parameter.
# HISTORY: Include trigger ID in subject of the "GRB not found in catalogs" email.
# HISTORY:
# HISTORY: Revision 1.14 2013/08/10 06:18:56 apsop
# HISTORY: Use FITSfile objects instead of system() calls and backticks to
# HISTORY: set/get the ANCRFILE and RESPFILE keywords (others already did).
# HISTORY: system() and backticks use the login shell environment, which is
# HISTORY: set up for an old headas and is incompatible with heasoft-6.13
# HISTORY: (eg. PFILES); FITSfile runs tools in the correct environment as set
# HISTORY: in sw0.par. Also added error handling on getting the RMF files, and
# HISTORY: list the xspec script. Added comments, esp. subroutine headers.
# HISTORY:
# HISTORY: Revision 1.13 2013/04/03 22:01:26 apsop
# HISTORY: sub sendEmail: strip apostrophes from subject, so they don't break
# HISTORY: the constructed command line.
# HISTORY:
# HISTORY: Revision 1.12 2012/08/30 07:36:30 apsop
# HISTORY: Many more log entries to tell which GRB coordinates are used and
# HISTORY: which catalog they came from, as in UvotProduct.pm and XrtProduct.pm.
# HISTORY: Reworded the email for Optical/XRT coordinates mismatch so it's clearer
# HISTORY: that is a warning and what should be checked. Note source of trigtime
# HISTORY: (although info is not currently used). Lots of minor formatting changes.
# HISTORY:
# HISTORY: Revision 1.11 2012/02/15 01:38:29 apsop
# HISTORY: The /* following setplot splashpage in the xspec command list should
# HISTORY: have been on its own line.
# HISTORY:
# HISTORY: Revision 1.10 2012/02/07 22:44:27 apsop
# HISTORY: Added object and sequence to "OPTICAL GRB coordinates
# HISTORY: outside XRT coordinates" alert message.
# HISTORY:
# HISTORY: Revision 1.9 2012/02/04 08:05:52 apsop
# HISTORY: Added "/*" to setplot splashpage off command, supposedly this avoids
# HISTORY: certain problems if used with an older version of Xspec.
# HISTORY:
# HISTORY: Revision 1.8 2012/02/02 06:22:39 apsop
# HISTORY: Added "setplot splashpage off" to xspec script.
# HISTORY:
# HISTORY: Revision 1.7 2012/01/12 06:52:03 apsop
# HISTORY: Changes going to proc3.15.03
# HISTORY:
# HISTORY: Revision 1.5 2011/01/20 19:00:12 apsop
# HISTORY: Added code to:
# HISTORY: 1 - use UVOT attitude file if available
# HISTORY: 2 - obtain correct TRIGTIME
# HISTORY: 3 - use correct burst RA and DEC
# HISTORY:
# HISTORY:
#
# VERSION: $Revision: 1.16 $
#
##############################################################################
use Subs::Sub;
use Util::GTIfile;
use Util::SWCatalogue;
use File::Path;
use File::Basename;
use Cwd;
@ISA = ("Subs::Sub");
use strict;
sub new {
my $proto=shift;
my $self=$proto->SUPER::new();
$self->{DESCRIPTION}="Make XRT GRB Light Curve";
return $self;
}
##################
# METHODS:
##################
sub body {
my $self=shift;
my $log =$self->log();
my $filename=$self->filename();
my $jobpar =$self->jobpar();
my $procpar = $self->procpar();
#########################
# make xrt grb light curve
#########################
$self->make_xrt_grb_lc();
} # end of body method
#############################################################################
#
# make_xrt_grb_lc: Produce xrt grb light curve files
#
# Parameters: none
#
#############################################################################
sub make_xrt_grb_lc {
my $self=shift;
my $log = $self->log();
my $filename = $self->filename();
my $jobpar = $self->jobpar();
my $procpar = $self->procpar();
my $watchers = $procpar->read("watchers");
my $object = $jobpar->read('object');
my $seq = $jobpar->read('sequence');
my $segnum = substr($seq, 8);
$segnum = $segnum*1;
if ($segnum > 0) {
$log->entry("Segment number is greater than zero: $segnum, no light curves or spectral analysis made");
return;
}
###################################
# CHECK IF WE HAVE ALL FILES NEEDED
###################################
# Event files for each mode. The PC and WT modes can be processed
# together, but the PD/LR and PD/PU modes must be done separately (and
# for each of B0 and B1). Note we only look for settled pointing (PO)
# files.
my @pcwpoFileList = $filename->get('event', 'xrt', 'pcw?po', '*');
my @wtwpoFileList = $filename->get('event', 'xrt', 'wtw?po', '*');
my @lrb0poFileList = $filename->get('event', 'xrt', 'lrb0po', '*');
my @lrb1poFileList = $filename->get('event', 'xrt', 'lrb1po', '*');
my @pub0poFileList = $filename->get('event', 'xrt', 'pub0po', '*');
my @pub1poFileList = $filename->get('event', 'xrt', 'pub1po', '*');
# get the output stem: sw or st and sequence number
my $outstem = undef;
foreach my $evt ( @wtwpoFileList, @pcwpoFileList,
@lrb0poFileList, @lrb1poFileList,
@pub0poFileList, @pub1poFileList ) {
$evt =~ /(s[wt]\d{11})x(pc|wt|lr|pu)/;
if ( $1 ) {
$outstem = $1;
last;
}
}
if (!defined $outstem) {
$log->error(1, "failed to determine outstem to run XrtGrbLc, exit 1");
return;
}
my $evtfiles = join(',', @pcwpoFileList, @wtwpoFileList); # all PC, WT modes
$log->entry("XRT PC & WT modes event files found: $evtfiles");
my $evtfiles_lrb0 = join(',', @lrb0poFileList); # PD/LR mode, B0
$log->entry("XRT PD/LR mode B0 event files found: $evtfiles_lrb0");
my $evtfiles_lrb1 = join(',', @lrb1poFileList); # PD/LR mode, B1
$log->entry("XRT PD/LR mode B1 event files found: $evtfiles_lrb1");
my $evtfiles_pub0 = join(',', @pub0poFileList); # PD/PU mode, B0
$log->entry("XRT PD/PU mode B0 event files found: $evtfiles_pub0");
my $evtfiles_pub1 = join(',', @pub1poFileList); # PD/PU mode, B1
$log->entry("XRT PD/PU mode B1 event files found: $evtfiles_pub1");
if ( (! $evtfiles) &&
(! $evtfiles_lrb0) && (! $evtfiles_lrb1) &&
(! $evtfiles_pub0) && (! $evtfiles_pub1) ) {
$log->error(1, "No event files exist, exit 1");
return ;
}
# xrthd, mkf files are used for PC, WT modes but not PD modes, so
# shouldn't return if they're not present as we used to.
my $xrthd_header = $filename->get('hk', 'xrt', 'hd');
my $have_xrthd = -e $xrthd_header;
if (! $have_xrthd ) {
$log->error(1, "$xrthd_header xrthd file does not exist, can't process XRT PC, WT modes");
}
my $mkf_filter = $filename->get('filter', 's');
my $have_mkf = -e $mkf_filter;
if (! $have_mkf ) {
$log->error(1, "$mkf_filter mkf file does not exist, can't process XRT PC, WT modes");
}
# Attitude file: pat if available, else sat.
# Return if none found, but it's not used for the PD modes so maybe we
# shouldn't?
my $attitude = $filename->get('attcorr', 'p');
if (-e $attitude) {
$log->entry("Got pat attitude file $attitude");
} else {
$log->error(1, "$attitude pat attitude file does not exist, trying sat file");
$attitude = $filename->get('attitude', 's');
if (-e $attitude) {
$log->entry("Got sat attitude file $attitude");
} else {
$log->error(1, "Unable to find attitude file, exit 1");
return ;
}
}
################################
# ADD KEYWORDS
################################
## (why do we have to go through all this to make the list of files??)
# my @evtfiles1 = split(/","/, "$evtfiles");
# my @evtfilesLR = split(/","/, "$evtfiles_lr");
# my @evtfilesPU = split(/","/, "$evtfiles_pu");
# my $evtfiles2 = join(',', @evtfiles1, @evtfilesLR, @evtfilesPU, $attitude);
#
# my @evtfiles3 = split(/,/, $evtfiles2);
#
# $self->write_standard_keyword(@evtfiles3);
$self->write_standard_keyword( ( @wtwpoFileList, @pcwpoFileList,
@lrb0poFileList, @lrb1poFileList,
@pub0poFileList, @pub1poFileList,
$attitude ) );
#################################ADD KEYWORDS
########################
# GET GRB COORDINATES
########################
# Read coords from job.par, but note these are never actually used:
# see comments below on choosing "best" coords.
my $ranom = $jobpar->read("burst_ra");
my $decnom = $jobpar->read("burst_dec");
# First try Lorella's catalog; if it fails,
# then try JD's catalog; then try local catalog;
# and finally send an email.
# Read local catalog
my $refPL = $self->Subs::UvotProduct::get_grb_coord_local();
# Check if this is a GRB, if not ($refPL == -1), simply return
# (Message is logged by get_grb_coord_local.)
if (defined $refPL and $refPL == -1){
return;
}
# Lorella's catalog.
my $refP = $self->Subs::UvotProduct::get_swiftgrb_coord();
my $catRead = "Lorella";
# JD's catalog
if (!defined $refP) {
$refP = $self->Subs::UvotProduct::get_grb_coord();
$catRead = "JD";
}
if (!defined $refP and !defined $refPL) {
my $enddate = $jobpar->read("enddate");
my $endtime = $jobpar->read("endtime");
my $dateobj = Util::Date->new($enddate, $endtime);
my $emjd = $dateobj->mjd();
my $ndateobj = Util::Date->new();
my $nmjd = $ndateobj->mjd();
my $alerT = $procpar->read("alerTime");
if (!defined $alerT ) {
$alerT = 2.0;
}
my $dt = $nmjd - $emjd;
if ($dt > $alerT) {
my $object = $jobpar->read('object');
my $trigid = $jobpar->read('sequence');
my $target = $jobpar->read('target');
my $subject = "GRB not found in catalogs: $trigid";
my $msg = "Unable to find GRB info in all catalogs for:\n" .
"object= $object\n" .
"sequence= $trigid\ntarget= $target\n\n" .
"\tPlease, make sure that such info is actually available.\n\n" .
"\tYou can add the required info to the local GRB catalog\n" .
"by using\n\n" .
"/data/sdc/local/data/sdc4/apsop/Processman/PMT/bin/AddGRB.pl\n\n";
$self->sendEmail($subject, $msg, $watchers);
}
$log->error(1, "Unable to find GRB in catalogs. Email notice sent to $watchers");
return;
}
my @rItems = qw/bat_ra bat_dec bat_pos_err xrt_ra xrt_dec xrt_pos_err
uvot_ra uvot_dec uvot_pos_err ot_ra ot_dec ot_pos_err/;
if (defined $refPL) {
if (!defined $refP) {
foreach my $it (@rItems) {
$refP->{$it} = $refPL->{$it};
}
$log->entry("Taking all coordinates from Local Catalog");
} else {
if ((defined $refPL->{bat_over} and $refPL->{bat_over} == 1) or
(!defined $refP->{bat_ra} or !defined $refP->{bat_dec})) {
$refP->{bat_ra} = $refPL->{bat_ra};
$refP->{bat_dec} = $refPL->{bat_dec};
$refP->{bat_pos_err} = $refPL->{bat_pos_err};
$log->entry("BAT invalid in $catRead cat, or Local override");
$log->entry("Using BAT coordinates from Local Catalog");
} else {
$log->entry("Using BAT coordinates from $catRead Catalog");
}
if ((defined $refPL->{xrt_over} and $refPL->{xrt_over} == 1) or
(!defined $refP->{xrt_ra} or !defined $refP->{xrt_dec})) {
$refP->{xrt_ra} = $refPL->{xrt_ra};
$refP->{xrt_dec} = $refPL->{xrt_dec};
$refP->{xrt_pos_err} = $refPL->{xrt_pos_err};
$log->entry("XRT invalid in $catRead cat, or Local override");
$log->entry("Using XRT coordinates from Local Catalog");
} else {
$log->entry("Using XRT coordinates from $catRead Catalog");
}
if ((defined $refPL->{ot_over} and $refPL->{ot_over} == 1) or
(!defined $refP->{ot_ra} or !defined $refP->{ot_dec})) {
$refP->{ot_ra} = $refPL->{ot_ra};
$refP->{ot_dec} = $refPL->{ot_dec};
$refP->{ot_pos_err} = $refPL->{ot_pos_err};
$log->entry("OPT invalid in $catRead cat, or Local override");
$log->entry("Using OPT coordinates from Local Catalog");
} else {
$log->entry("Using OPT coordinates from $catRead Catalog");
}
if ((defined $refPL->{uvot_over} and $refPL->{uvot_over} == 1) or
(!defined $refP->{uvot_ra} or !defined $refP->{uvot_dec})) {
$refP->{uvot_ra} = $refPL->{uvot_ra};
$refP->{uvot_dec} = $refPL->{uvot_dec};
$refP->{uvot_pos_err} = $refPL->{uvot_pos_err};
$log->entry("UVOT invalid in $catRead cat, or Local override");
$log->entry("Using UVOT coordinates from Local Catalog");
} else {
$log->entry("Using UVOT coordinates from $catRead Catalog");
}
}
} else {
$log->entry("Not in Local Catalog. Taking all coordinates from $catRead Catalog");
}
# The GRB position must be chosen based on the "best" coordinates. If UVOT
# coordinate are defined they will be used, else the Optical or OT
# coordinates will be used if defined (with a warning email sent if they
# disagree too much with the XRT ones), and finally if the XRT coordinates
# are defined they will be used. No light curve is made if only BAT
# coordinates are available.
if (defined $refP->{uvot_ra} && defined $refP->{uvot_dec} &&
(!defined $refPL or !defined $refPL->{uvot_bad_data} or
$refPL->{uvot_bad_data} == 0)) {
$ranom = $refP->{uvot_ra};
$decnom = $refP->{uvot_dec};
$log->entry("Will be using UVOT coordinates ($ranom, $decnom)");
} elsif (defined $refP->{ot_ra} && defined $refP->{ot_dec} &&
(!defined $refPL or !defined $refPL->{ot_bad_data} or
$refPL->{ot_bad_data} == 0)) {
$ranom = $refP->{ot_ra};
$decnom = $refP->{ot_dec};
$log->entry("Will be using OPTICAL coordinates ($ranom, $decnom)");
if (defined $refP->{xrt_ra} && defined $refP->{xrt_dec} &&
defined $refP->{xrt_pos_err}) {
my $err = 2.0*$refP->{xrt_pos_err}/3600.0;
my $dra = $refP->{xrt_ra}-$err;
my $ura = $refP->{xrt_ra}+$err;
my $ddec = $refP->{xrt_dec}-$err;
my $udec = $refP->{xrt_dec}+$err;
if (($ranom < $dra or $ranom > $ura) or
($decnom < $ddec or $decnom > $udec)) {
my $subject = "OPTICAL GRB coordinates outside XRT coordinates";
my $msg = "Object = $object\nSequence = $seq\n\n"
. "Warning: While using OPTICAL coordinates for GRB, it appears that\n"
. "these coordinates ($ranom, $decnom) are not in agreement\n"
. "with the XRT ones ("
. $refP->{xrt_ra} . ',' . $refP->{xrt_dec} . ") within 2.*xrt_pos_err of "
. $refP->{xrt_pos_err} . " arcsec.\n"
. "\n\tPlease, verify that the catalog values are as expected.\n"
. "\tPositions from $catRead catalog unless Local overrides.\n";
$self->sendEmail($subject, $msg, $watchers);
$log->entry("make_xrt_grb_lc: Warning sent for Optical/XRT disagreement > 2*xrt_pos_err\n"
. " xrt_ra=" . $refP->{xrt_ra} . " deg, xrt_dec=" . $refP->{xrt_dec}
. " deg, xrt_pos_err=" . $refP->{xrt_pos_err} . " arcsec" );
}
}
} elsif (defined $refP->{xrt_ra} && defined $refP->{xrt_dec} &&
(!defined $refPL or !defined $refPL->{xrt_bad_data} or
$refPL->{xrt_bad_data} == 0)) {
$ranom = $refP->{xrt_ra};
$decnom = $refP->{xrt_dec};
$log->entry("Will be using XRT coordinates ($ranom, $decnom)");
} else {
$log->error(1, "Only BAT coordinates, or none, are available. No XRT light curve for GRB will be made.");
return;
}
my $plotfile0 = $filename->get("lcplot", "xrt", "sr", 0);
my $plotfile1 = $filename->get("lcplot", "xrt", "sr", 1);
my $plotfile2 = $filename->get("lcplot", "xrt", "sr", 2);
#############################
#
# MAKE XRT GRB LIGHT CURVE
# PC AND WT MODES
#
#############################
if ( $evtfiles && $have_xrthd && $have_mkf ) {
$log->entry('Running xrtgrblc for PC & WT modes');
my $odir = './tmp_dir_pcwt';
if (-e $odir) {
system("rm -fr $odir");
}
# xrtgrblc needs one of these for every evtfile, even if they're all
# the same.
my $num1 = $#pcwpoFileList + 1;
my $num2 = $#wtwpoFileList + 1;
my $mkffiles = join(',', ($mkf_filter)x($num1 + $num2));
my $xhdfiles = join(',', ($xrthd_header)x($num1 + $num2));
my $attfiles = join(',', ($attitude)x($num1 + $num2));
my $xrtgrblc=Util::HEAdas->new("xrtgrblc")->is_script( 1 );
$xrtgrblc->params({
evtfiles => $evtfiles,
mkffiles => $mkffiles,
xhdfiles => $xhdfiles,
attfiles => $attfiles,
outdir => $odir,
outstem => $outstem,
srcra => $ranom,
srcdec => $decnom,
minenergy => 0.3,
maxenergy => 10.0,
splitenergy => 0.0,
splitorbits => 'no',
cutlowbins => 'yes',
plotftype => '/gif',
clobber => 'yes',
history => 'yes',
cleanspec => 'yes'
});
$xrtgrblc->run();
# Replot the light curve. A better title should be added, and
# the X-axis label should not say "BAT Trigger" if no BAT trigger
# Setup some plot labels
my $trigtime = $jobpar->read('burst_time');
my $trigtimeSrc = "burst_time Parameter";
my $trigdate = undef;
my $xlab;
# Note: trigger time read here appears to only be used for the axis label
# of the xsr_lc plot below!
if (!defined $trigtime or $trigtime == 0){
$trigtime = $self->Subs::UvotProduct::getTrigFromDB();
$trigtimeSrc = "TDRSS Catalog";
if (!defined $trigtime){
$trigtime = $self->Subs::UvotProduct::getTrigFromSWIFTCatalog();
$trigtimeSrc = "Lorella Catalog";
if (!defined $trigtime){
$trigtime = $self->Subs::UvotProduct::getTrigFromJDCatalog();
$trigtimeSrc = "JD Catalog";
if (!defined $trigtime){
$trigtime = $self->Subs::UvotProduct::getTrigFromLocalCatalog();
$trigtimeSrc = "Local Catalog";
}
}
}
}
if (defined $trigtime){
my $dateobj = Util::Date->new($trigtime);
$trigdate = $dateobj->date().'T'.$dateobj->time();
}
if (defined $trigtime and defined $trigdate){
$xlab = "Time(s) since BAT trigger time(UT $trigdate/MET $trigtime)";
} else {
$xlab = "Time [MET(s)]";
}
my $evtfile = ( @pcwpoFileList, @wtwpoFileList )[ 0 ];
my $evtfits = Util::FITSfile->new($evtfile);
my $dateobs = $evtfits->keyword( 'DATE-OBS' );
$dateobs =~ s/[T ]\d{2}:\d{2}:[\d\.]+$//;
my $title = uc( $object ) . ' SWIFT XRT '
. $self->jobpar->read('sequence')
. ' (' . $dateobs . ')';
# change to the dir where the qdp file is
my $curdir = getcwd( );
chdir $odir;
my $qdp = Util::Xanadu->new( 'qdp' );
$qdp->verbose( 1 );
$qdp->seriousness( 1 );
$qdp->script(
"${outstem}_xpwetsrab.qdp",
'/null',
'SCR WHITE',
'TIME OFF',
'LAB F',
"LAB T \"$title\"",
"LAB X \"$xlab\"",
'YPL OFF 1 3 4 5 7 8',
'LOG',
'GAP ER',
'RES',
"cpd ${outstem}xsr_lc.gif/gif",
"PLOT",
"QUIT"
)->run( );
chdir $curdir;
#
# plot the spectra
#
my @phaF = ();
push @phaF, "$odir/${outstem}_xpcetsr.pha";
push @phaF, "$odir/${outstem}_xwtetsr.pha";
$self->plotSpectra( \@phaF, "${outstem}xsr_ph.gif/gif" ); # NOTE basename here
#
# move files into their final names
# sr = source
#
rename "$odir/${outstem}_xpcetsr.pha", "$odir/${outstem}xpcsr.pha";
rename "$odir/${outstem}_xpcet.arf", "$odir/${outstem}xpcsr.arf";
rename "$odir/${outstem}_xwtetsr.pha", "$odir/${outstem}xwtsr.pha";
rename "$odir/${outstem}_xwtet.arf", "$odir/${outstem}xwtsr.arf";
rename "$odir/${outstem}_xpcetsra.lc", "$odir/${outstem}xpcsr.lc";
rename "$odir/${outstem}_xwtetsra.lc", "$odir/${outstem}xwtsr.lc";
# remove un-needed files
unlink "$odir/${outstem}_xpcetbg.pha", "$odir/${outstem}_xwtetbg.pha",
"$odir/${outstem}_xpwetsrfb_lc.gif", "$odir/${outstem}_xpwetsrcb_lc.gif",
"$odir/${outstem}_info.fits";
#unlink "$odir/${outstem}_xpwetsrab.qdp", "$odir/${outstem}_xpwetsrab.pco";
# Fixup keywords: rewrite (or add) ANCRFILE with the arf file's new name,
# on both HDUs.
if ( -e "$odir/${outstem}xpcsr.pha" ) {
Util::FITSfile->new("$odir/${outstem}xpcsr.pha", 0)
->keyword("ANCRFILE", "${outstem}xpcsr.arf");
Util::FITSfile->new("$odir/${outstem}xpcsr.pha", 1)
->keyword("ANCRFILE", "${outstem}xpcsr.arf");
}
if ( -e "$odir/${outstem}xwtsr.pha" ) {
Util::FITSfile->new("$odir/${outstem}xwtsr.pha", 0)
->keyword("ANCRFILE", "${outstem}xwtsr.arf");
Util::FITSfile->new("$odir/${outstem}xwtsr.pha", 1)
->keyword("ANCRFILE", "${outstem}xwtsr.arf");
}
# Concat log files with the ones in the main directory
$self->concatLogs( $outstem, $odir, $curdir );
# Move any remaining files in $dir to $curdir.
if (-e $odir) {
system("mv $odir/* $curdir");
system("rm -fr $odir");
}
} # end of PC, WT modes
##########################################################################
#
# MAKE XRT GRB LIGHT CURVES
# PD MODES (LR and PU)
#
# Names for the output files try to follow the PC,WT convention above
# (eg, "${outstem}xpcsr.lc") but need to include b0/b1. PC,WT names
# include "sr" ("Swift GRB Product File Naming Convention" says
# sr=source); even though xrtproducts help says PD modes have no spatial
# information and all events are used, Matteo Perri (XRT team at ASDC)
# says "sr" is appropriate for these modes too because they were only
# used for extremely bright sources. "po" is not included because we're
# using only the settled pointing event files.
##########################################################################
if ( $evtfiles_lrb0 ) {
$self->make_pd_mode_lc( $evtfiles_lrb0, 'xlrb0sr', $outstem );
}
if ( $evtfiles_lrb1 ) {
$self->make_pd_mode_lc( $evtfiles_lrb1, 'xlrb1sr', $outstem );
}
if ( $evtfiles_pub0 ) {
$self->make_pd_mode_lc( $evtfiles_pub0, 'xpub0sr', $outstem );
}
if ( $evtfiles_pub1 ) {
$self->make_pd_mode_lc( $evtfiles_pub1, 'xpub1sr', $outstem );
}
} # end of make_xrt_grb_lc method
##########################################################################
#
# make_pd_mode_lc( $evtfiles, $basename, $outstem )
#
# Run xrtproducts to make light curve, arf, and pha products for one of the
# XRT Photo Diode (PD/LR and PD/PU) modes: xlrb0, xlrb1, xpub0,
# xpub1. (b1=bias subtracted onboard, b0=not subtracted.) They're all done
# the same way, but unlike PC and WT must be done separately; call this
# with the appropriate event file for each one. Files created are *.lc,
# *.arf, *.pha, *_lc.gif, and *_ph.gif. Work is done in a temporary
# directory, as for the PC, WT modes. The procedure used is from Matteo
# Perri (ASDC) (emails 2014-03-28 and 2014-04-07).
#
# Parameters:
# $evtfiles: String of input event file names for this mode, comma
# separated if more than one.
# $basename: Base name of the output files, will be appended onto $outstem.
# Eg., 'xlrb0sr'.
# $outstem: Stem for the output file names. Eg, 'sw00106415000'.
#
# Eg., $basename='xlrb0sr', $outstem='sw00106415000' produces files named
# sw00106415000xlrb0sr.lc, etc.
#
##########################################################################
sub make_pd_mode_lc {
my $self = shift;
my ($evtfiles, $basename, $outstem) = @_;
my $log = $self->log();
$log->entry("Running xrtproducts for PD mode ${basename}");
my $curdir = getcwd( );
my $odir = "./tmp_dir_${basename}";
if (-e $odir) {
system("rm -fr $odir");
}
# Names for the output files.
my $lcName = "${outstem}${basename}.lc";
my $phaName = "${outstem}${basename}.pha";
my $arfName = "${outstem}${basename}.arf";
# Run xrtproducts to create the light curve and related files. The
# parameters listed were specified by Matteo Perri (mostly) as the
# appropriate values for processing all the PD modes and should be the
# complete set (email from Matteo 2014-04-07).
# NOTE: Setting stemout=>"${outstem}${basename}" is a partial
# workaround suggested by Matteo for a bug in xrtproducts v.0.4.1 and
# earlier, which hardcodes the plot title in the light curve gif to
# "Light curve (${stemout}pcsr.lc)" for ALL modes. (The pha gif title
# appears to use $phafile.) This also changes the names of the .gif
# files, but since we specify all the others explicitly, it appears to
# have no other effect. Change this back to stemout=>$outstem, and the
# rename commands below, if Matteo gives us a patch for the bug.
my $xrtproducts = Util::HEAdas->new("xrtproducts")->is_script( 1 );
$xrtproducts->params( {
infile => $evtfiles,
outdir => $odir,
lcfile => $lcName,
phafile => $phaName,
arffile => $arfName,
# This stemout partially works around the bug in lc.gif plot titles
stemout => "${outstem}${basename}",
binsize => -99,
display => 'no',
gtifile => 'NONE',
rmffile => 'CALDB',
mirfile => 'CALDB',
transmfile => 'CALDB',
inarfile => 'CALDB',
psfile => 'CALDB',
vigfile => 'CALDB',
psfflag => 'yes',
extended => 'no',
plotdevice => 'gif',
pilow => 30,
pihigh => 1000,
clobber => 'yes',
chatter => 3,
history => 'yes',
cleanup => 'yes'
} );
$xrtproducts->run();
# Move files into their final names
# Doesn't look like xrtproducts has a way to set the .gif names;
# xrtproducts names the gif files ${stemout}lc.gif and ${stemout}ph.gif
rename "${odir}/${outstem}${basename}lc.gif", "${odir}/${outstem}${basename}_lc.gif";
rename "${odir}/${outstem}${basename}ph.gif", "${odir}/${outstem}${basename}_ph.gif";
# remove unneeded files
# Fixup ANCRFILE keywords in the .pha file with the .arf file's name
if ( -e "${odir}/${phaName}" ) {
my $phafile = Util::FITSfile->new("${odir}/${phaName}");
$phafile->ext(0);
$phafile->keyword("ANCRFILE", $arfName);
$phafile->ext(1);
$phafile->keyword("ANCRFILE", $arfName);
}
# Concat log files with the ones in the main directory
$self->concatLogs( $outstem, $odir, $curdir );
# Move any remaining files in $dir to $curdir.
if (-e $odir) {
system("mv $odir/* $curdir");
system("rm -fr $odir");
}
} # end of make_pd_mode_lc method
###############################################################
#
# write_standard_keyword: Add standard keywords to FITS files.
# This one is currently called by XrtGrbLc and UvotProducts.
# There is a very similar routine in SW0WrapUp (this one was probably
# copied from there).
#
# Parameter:
# @file1: List of filenames to add keywords to.
#
###############################################################
sub write_standard_keyword {
my $self=shift;
my $log = $self->log();
my $filename = $self->filename();
my $jobpar = $self->jobpar();
my $procpar = $self->procpar();
my @file1 = @_;
#######################
# Setup for UTCF value
#######################
my $time_file = $filename->get('timedata', 'swift', '', 0);
my (@utcf_times, @utcf);
if( -f $time_file ){
my $time_fits = Util::FITSfile->new($time_file, 'UTCF');
@utcf_times = $time_fits->cols('TIME')->table();
@utcf = $time_fits->cols('UTCF')->table();
}
my $soft_version = $jobpar->read('softver');
my $cal_version = $jobpar->read('caldbver');
my $reprocess = $jobpar->read('reprocess');
############################
# Trigger time, if relevant
############################
my $trigtime = $jobpar->read('burst_time');
####################################
# loop over all FITS files
####################################
my $file;
foreach $file (@file1 ) {
# $log->entry("Doing file $file.");
$log->entry("Adding standard keywords to file $file");
my $fits = Util::FITSfile->new($file);
my $is_tdrss = $file =~ /s[wt][t\d]\d{10}ms/;
my $is_batevt = $file =~ /s[wt][t\d]\d{10}bev/;
my $is_eng = $file =~ /s[wt][t\d]\d{10}.*en\.hk$/;
my $is_att = $file =~ /s[wt][t\d]\d{10}.at\.fits$/;
my $is_sm = $file =~ /s[wt][t\d]\d{10}bsmcb\.fits$/;
##################################
# loop over HDUs in the FITS file
##################################
my $nhdus = $fits->nhdus();
unless ($nhdus) {
# $log->error(2, "Cannot get number of HDUs for FITS file $file.");
next;
}
my $hdu;
for($hdu=0; $hdu<$nhdus; $hdu++) {
$fits->ext($hdu);
my $extname = $hdu==0 ? '' : $fits->keyword('EXTNAME');
################################
# write keywords to the file
################################
$fits->begin_many_keywords();
$fits->keyword('PROCVER', $procpar->read('version'),
'Processing script version' );
unless( $is_eng || ($extname && $extname=~/GTI/) ){
$fits->keyword('SOFTVER', $soft_version, 'HEASOFT and Swift versions');
$fits->keyword('CALDBVER', $cal_version, 'CALDB index versions used');
}
if($hdu==0){
$fits->keyword('TIMESYS', 'TT', 'time system');
$fits->keyword('MJDREFI', 51910, 'MJD reference day 01 Jan 2001 00:00:00');
$fits->keyword('MJDREFF', 7.428703700000000E-04,
'MJD reference (fraction of day) 01 Jan 2001 00:00:00');
$fits->keyword('CLOCKAPP', 'F', 'If clock correction are applied (F/T)');
$fits->keyword('TIMEUNIT', 's', 'Time unit for timing header keywords');
$fits->keyword('REPROC', 'T', 'Is this from a bulk reprocessing run?') if $reprocess eq 'yes';
}else{
$fits->keyword('TIERRELA', '1.0E-8', '[s/s] relative errors expressed as rate');
$fits->keyword('TIERABSO', '1.0', '[s] timing precision in seconds');
}
unless($is_tdrss || $is_sm){
$fits->keyword('OBS_ID', sprintf("'%011d'", $jobpar->read('sequence')),
'Observation ID' );
$fits->keyword('SEQPNUM', int($jobpar->read('seqprocnum')),
'Number of times the dataset processed' );
$fits->keyword('TARG_ID', int(sprintf("%d",$jobpar->read('target'))),
'Target ID');
$fits->keyword('SEG_NUM', int(sprintf("%d",$jobpar->read('obs'))),
'Segment number');
$fits->keyword('OBJECT', $jobpar->read('object'), 'Object name');
$fits->keyword('RA_OBJ', $jobpar->read('burst_ra'), '[deg] R.A. Object');
$fits->keyword('DEC_OBJ', $jobpar->read('burst_dec'), '[deg] Dec Object');
$fits->keyword('RA_PNT', $jobpar->read('ra'), '[deg] RA pointing');
$fits->keyword('DEC_PNT', $jobpar->read('dec'), '[deg] Dec pointing');
$fits->keyword('PA_PNT', $jobpar->read('roll'), '[deg] Position angle (roll)');
$fits->keyword('TRIGTIME', $trigtime, '[s] MET TRIGger Time for Automatic Target')
if $trigtime;
}
if($is_batevt){
$fits->keyword('CATSRC',
$jobpar->read('burst_cat_src') eq 'yes' ? 'T' : 'F');
}
if ( $is_att ) {
$fits->keyword('-UTCFINIT', ' ');
$fits->keyword('-ATTSTATU', ' ');
} else {
#############################
# Determine and set UTCFINIT
#############################
my $tstart = $fits->keyword('TSTART');
unless( $tstart ){
my $date = $fits->keyword('DATE-OBS');
if( $date ){
my $start = Util::Date->new($date);
$tstart = $start->seconds();
}
}
if( $tstart && @utcf ){
my $itime = 0;
while( $tstart > $utcf_times[$itime] &&
$itime < $#utcf_times){ $itime++ }
$fits->keyword('UTCFINIT', $utcf[$itime], '[s] UTCF at TSTART');
}
}
$fits->end_many_keywords();
}
}
} # end of write_standard_keyword method
###############################################################
#
# sendEmail: Send an email
#
# Parameters:
# $subject: Subject
# $msg: Text of the message, can contain \n etc.
# $watchers: Recipient(s), email addresses separated by spaces
#
# Apostrophes will be removed from subject, because they break the
# shell command line we construct. They're OK within the msg text.
###############################################################
sub sendEmail {
my $self=shift;
my ($subject, $msg, $watchers) = @_;
my $jobpar = $self->jobpar();
my $trigid =$jobpar->read('sequence');
my $tm = time();
my $tfile = "/tmp/msg_$tm".'_'.$trigid;
# Remove any apostrophes from subject, because
# they will break the constructed command line.
$subject =~ s/\'//g ;
if (open OUTF, ">$tfile") {
print OUTF "$msg";
close OUTF;
my $hostname = `hostname`;
chomp $hostname;
if ($hostname =~ /^sdc/ ) {
my $cmd = "mail -s \'$subject\' $watchers < $tfile";
`$cmd`;
unlink $tfile;
} else {
my $rv = system("scp -q $tfile apsop\@sdc:/tmp >& /dev/null");
if ($rv == 0) {
my $cmd1 = "ssh -nq apsop\@sdc \"mail -s \'$subject\' $watchers < $tfile\"";
my $rval = system($cmd1);
unlink $tfile;
system("ssh -nq apsop\@sdc rm -f $tfile");
}
}
}
} # end of sendEmail method
###########################################################
#
# plotSpectra - groups and plots spectra.
# Runs grppha and xspec.
#
# Parameters:
# $refspectra: Reference to array of .pha files to plot.
# $outplot: Name of the resulting plot file. Deleted if
# no plot is made (ie, plot file comes out empty).
#
############################################################
sub plotSpectra {
my $self = shift;
my ( $refspectra, $outplot ) = @_;
my $log = $self->log();
my $jobpar = $self->jobpar();
my $procpar = $self->procpar();
my $watchers = $procpar->read("watchers");
my $headas = $procpar->read('headas');
my $GRB = $jobpar->read('object');
$GRB = uc($GRB);
my @localclean = ( ); # list of files we'll delete before returning
# check which spectra we actually have
my @spectra = ( );
foreach my $spec1 ( @$refspectra ) {
if ( -e $spec1 ) {
push @spectra, basename( $spec1 );
}
}
return unless @spectra;
# change to the dir where the spectra are
my $workdir = dirname( $refspectra->[ 0 ] );
my $curdir = getcwd( );
chdir $workdir;
#
# Get the response matrices and group each spectrum
#
my $i = 1;
my @grpspectra = ( );
my $grppha=Util::HEAdas->new("grppha");
my @color = ();
foreach my $spec ( @spectra ) {
# Get name of response matrix file (RMF) from RESPFILE keyword
# (FITSfile->keyword strips quotes and whitespace (inc. \n))
# If can't, throw error and go on to next spectrum.
my $resp = Util::FITSfile->new($spec,1)->keyword('RESPFILE');
if ( $resp ) { # defined and not null
$log->entry("XrtGrbLc plotSpectra: $spec RESPFILE is $resp\n");
} else {
$log->error(2, "Could not read keyword RESPFILE from $spec");
next;
}
# Copy RMF file from caldb.
# Error and go on to next if we can't (there may be more info in the
# proclog file).
# Should be more general than this, but RMF should be in current
# directory.
my $caldbResp = "$ENV{CALDB}/data/swift/xrt/cpf/rmf/$resp";
if ( -e $caldbResp ) {
my $trv = system( "cp $caldbResp ./$resp" );
if ($trv == 0) {
$log->entry("got RMF file $resp from CALDB.\n");
push @localclean, $resp;
} else {
$log->error(2,"Could not copy $caldbResp to $workdir, status=$trv");
next;
}
} else {
$log->error(2,"RMF file $caldbResp not found in CALDB!");
next;
}
# Run grppha to group
my $grpspec = $spec;
$grpspec =~ s/\.(pha|pi)$/_grp.$1/;
$grppha->params( {
infile => $spec,
outfile => $grpspec,
comm => 'group min 20',
tempc => 'exit',
clobber => 'yes'
} );
$grppha->run();
push @grpspectra, $grpspec;
push @localclean, $grpspec;
# Build up color command for plotting
$grpspec =~ /x(pc|wt)/;
my $mode = $1;
push @color, sprintf( "setplot command color %d on %d",
$mode eq 'pc' ? 2 : 4, $i );
$i++;
}
#
# Write and run the xspec script
# A prettier plot should be created (using setplot command)
#
# Setup some plot labels
my $fits = Util::FITSfile->new($grpspectra[ 0 ]);
my $dateobs = $fits->keyword( 'DATE-OBS' );
$dateobs =~ s/[T ]\d{2}:\d{2}:[\d\.]+$//;
my $title = $GRB . ' SWIFT XRT '
. $self->jobpar->read('sequence')
. ' (' . $dateobs . ')';
my $xspec = Util::Xanadu->new("xspec");
$xspec->verbose( 1 );
$xspec->seriousness( 1 );
# Build the xspec command list
my @commands = ();
push @commands, ('query no',
'setplot splashpage off',
'/*',
'setplot energy',
'setplot command LA OT " "',
"setplot command LA T $title",
'setplot command TIME OFF',
'setplot command SCR WHITE',
'setplot command GAP ER',
'setplot command RES',
);
foreach my $c (@color) {
if ($c !~ /^\s*$/) {
push @commands, $c;
}
}
$i = 1;
foreach my $spec ( @grpspectra ) {
push @commands, ("data $i:$i $spec",
'/*',
'/*',
'ignore bad',
"ignore $i:**-0.3 10.0-**",
);
$i++;
}
push @commands, ("cpd $outplot",
'plot ldata',
'exit'
);
$log->entry("XrtGrbLc plotSpectra: xspec command list:\n");
$log->text( join("\n", @commands) );
# Execute xspec
$xspec->script(@commands)->run();
my $gif = (split/\//, $outplot)[0];
if (!-s $gif) {
$log->error(1, "Attempted to create plot file $gif but its size is 0. File removed");
unlink $gif;
}
# Remove files no longer needed and return to main directory
unlink @localclean;
chdir $curdir;
} # end of plotSpectra method
##################################################################
#
# concatLogs( $outstem, $odir, $curdir )
#
# Concatenate the log files in $odir (./tmp_dir_*) to the ones in the main
# working directory $curdir. (Output buffering should not be a problem,
# because Log.pm opens and closes the files for each entry written.) The
# file names all begin with $outstem. Removes the log files from $odir.
#
################################################################
sub concatLogs {
my $self = shift;
my ($outstem, $odir, $curdir) = @_;
if (-e $odir) {
if (-e "$odir/${outstem}pjl.html"){
`cat $odir/${outstem}pjl.html >> $curdir/${outstem}pjl.html`;
unlink "$odir/${outstem}pjl.html";
}
if (-e "$odir/${outstem}per.html"){
`cat $odir/${outstem}per.html >> $curdir/${outstem}per.html`;
unlink "$odir/${outstem}per.html";
}
if (-e "$odir/${outstem}pin.html"){
`cat $odir/${outstem}pin.html >> $curdir/${outstem}pin.html`;
unlink "$odir/${outstem}pin.html";
}
}
return;
} # end of concatLogs method
1