package Subs::XrtGrbLc;
##############################################################################
#
# DESCRIPTION: Produce XRT GRB lightcurve gif file
#
# HISTORY:
# HISTORY: $Log: XrtGrbLc.pm,v $
# 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: 0.0
#
##############################################################################
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 gif 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##############
my @pcwpoFileList = $filename->get('event', 'xrt', 'pcw?po','*');
my @wtwpoFileList = $filename->get('event', 'xrt', 'wtw?po','*');
# get the output stem
my $outstem = undef;
foreach my $evt ( @wtwpoFileList, @pcwpoFileList ) {
$evt =~ /(s[wt]\d{11})x(pc|wt)/;
if ( $1 ) {
$outstem = $1;
last;
}
}
if (!defined $outstem) {
$log->error(1, "failed to determine outstem to run XrtGrbLc, exit 1");
}
my $evtfiles = join(',', @pcwpoFileList, @wtwpoFileList);
if (! $evtfiles ) {
$log->error(1, "$evtfiles evt file not exist, exit 1");
return ;
}
my $xrthd_header = $filename->get('hk', 'xrt', 'hd');
if (! -e $xrthd_header ) {
$log->error(1, "$xrthd_header xrthd file not exist, exit 1");
return ;
}
my $mkf_filter = $filename->get('filter', 's');
if (! -e $mkf_filter ) {
$log->error(1, "$mkf_filter mkf file not exist, exit 1");
return ;
}
my $sat_attitude = $filename->get('attcorr', 'u');
if (! -e $sat_attitude) {
$log->error(1, "$sat_attitude UVOT attitude file not exist, trying pat file");
$sat_attitude = $filename->get('attcorr', 'p');
if(! -e $sat_attitude) {
$sat_attitude = $filename->get('attitude', 's');
if (! -e $sat_attitude) {
$log->error(1, "$sat_attitude sat file not exist, exit 1");
return ;
}
}
}
#################################ADD KEYWORDS
my @evtfiles1=split(/","/,"$evtfiles");
my $evtfiles2 = join(',', @evtfiles1, $sat_attitude);
my @evtfiles3=split(/,/,$evtfiles2);
$self->write_standard_keyword(@evtfiles3);
#################################ADD KEYWORDS
# Read coords from job.par, but note these are never actually used:
# see section below on choosing "best" coords.
my $ranom = $jobpar->read("burst_ra");
my $decnom = $jobpar->read("burst_dec");
##########GET GRB COORD
# First try Lorella catalog, if it fails,
# then try JD 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";
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 $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(',', ($sat_attitude)x($num1 + $num2));
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##############
my $odir = './tmp_dir';
if (-e $odir) {
system("rm -fr $odir");
}
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
#
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");
}
#
# Concatenate the log files in $odir (./tmp_dir) to the ones in the
# main working directory. (Output buffering should not be a problem,
# because Log.pm opens and closes the files for each entry written.)
# Then move any remaining files in ./tmp_dir to the main directory.
#
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";
}
system("mv $odir/* $curdir");
system("rm -fr $odir");
}
}# end of make_xrt_grb_lc method
###############################################################
#
# write_standard_keyword: Add standard keywords to all FITS files
#
# 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 = @_;
# $log->entry("Writing standard keywords to all FITS files");
#######################
# 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;
my $sat_attitude = $filename->get('attcorr', 'u');
if (! -e $sat_attitude) {
$sat_attitude = $filename->get('attcorr', 'p');
if (! -e $sat_attitude) {
$sat_attitude = $filename->get('attitude', 's');
}
}
my $attflag = '100';
if($sat_attitude =~ /uat\.fits/){
$attflag = '111';
} elsif ($sat_attitude =~ /pat\.fits/){
$attflag = '110';
}
foreach $file (@file1 ) {
# $log->entry("Doing 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$/;
####################################
# Which attitude file is being used
####################################
##################################
# 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';
# my $att0 = $fits->keyword('ATTFLAG');
# $attflag = $att0 if $att0;
}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{
$fits->keyword('ATTFLAG', "'$attflag'", 'Orgin of attitude information');
#############################
# 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");
}
}
}
}
###########################################################
#
# 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;
}
1