package Subs::XrtGrbLc;
##############################################################################
#
# DESCRIPTION: Produce XRT GRB lightcurve gif file
#
# HISTORY:
# HISTORY: $Log: XrtGrbLc.pm,v $
# 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
#############################################################################
# produce xrt grb light curve gif files
#############################################################################
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
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.
my $refPL = $self->Subs::UvotProduct::get_grb_coord_local(); # read local catalog
if(defined $refPL and $refPL == -1){
return;
}
my $refP = $self->Subs::UvotProduct::get_swiftgrb_coord(); # Lorella catalog
if(defined $refP and $refP == -1){
return;
}
if (!defined $refP) {
$refP = $self->Subs::UvotProduct::get_grb_coord(); # JD catalog
}
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:\nobject= $object\nsequence= $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\nby using\n\n/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};
}
} 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("Using BAT coordinates from Local 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("Using XRT coordinates from Local 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("Using OPT coordinates from Local 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("Using UVOT coordinates from Local Catalog");
}
}
}
# The GRB position must be chosen based on the "best" coordinates
# If UVOT coordinate are defined they will be used, the the Optical
# or OT coordinates will be used if defined, and finally if the
# XRT coordinates are defined they will be used.
# No light curve is made if only BAT corinates 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 = 'While using OPTICAL coordinates for GRB, it appears that'."\nthese coordinates ($ranom, $decnom) are not in agreement\nwith the XRT ones (".$refP->{xrt_ra}.','.$refP->{xrt_dec}.") with pos_err of (".$refP->{xrt_pos_err}.") arsec.\n\n\tPlease, verify.";
$self->sendEmail($subject, $msg, $watchers);
}
}
} 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 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 $trigdate = undef;
my $xlab;
if(!defined $trigtime or $trigtime == 0){
$trigtime = $self->Subs::UvotProduct::getTrigFromDB();
if(!defined $trigtime){
$trigtime = $self->Subs::UvotProduct::getTrigFromSWIFTCatalog();
if(!defined $trigtime){
$trigtime = $self->Subs::UvotProduct::getTrigFromJDCatalog();
if(!defined $trigtime){
$trigtime = $self->Subs::UvotProduct::getTrigFromLocalCatalog();
}
}
}
}
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
if ( -e "$odir/${outstem}xpcsr.pha" ) {
system( "fparkey value=\"${outstem}xpcsr.arf\" "
. "fitsfile=\"$odir/${outstem}xpcsr.pha\[0\]\" "
. "keyword=\"ANCRFILE\" add=yes" );
system( "fparkey value=\"${outstem}xpcsr.arf\" "
. "fitsfile=\"$odir/${outstem}xpcsr.pha\[1\]\" "
. "keyword=\"ANCRFILE\" add=yes" );
}
if ( -e "$odir/${outstem}xwtsr.pha" ) {
system( "fparkey value=\"${outstem}xwtsr.arf\" "
. "fitsfile=\"$odir/${outstem}xwtsr.pha\[0\]\" "
. "keyword=\"ANCRFILE\" add=yes" );
system( "fparkey value=\"${outstem}xwtsr.arf\" "
. "fitsfile=\"$odir/${outstem}xwtsr.pha\[1\]\" "
. "keyword=\"ANCRFILE\" add=yes" );
}
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
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
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;
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
#
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 = ( );
# 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 ) {
my $resp = `fkeypar "$spec+1" RESPFILE ; pget fkeypar value`;
chomp $resp;
$resp =~ s/['"]//g;
# should be more general than this, but RMF should be in current
# directory
if ( -e "$ENV{CALDB}/data/swift/xrt/cpf/rmf/$resp" ) {
# symlink "$ENV{CALDB}/data/swift/xrt/cpf/rmf/$resp", $resp;
my $trv = system( "cp $ENV{CALDB}/data/swift/xrt/cpf/rmf/$resp ./$resp");
if ($trv == 0) {
push @localclean, $resp;
}
}
# 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 );
my @commands = ();
push @commands, ('query no',
'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'
);
# my $script = "tmp.xspec.$$.xsl";
# my $str = join("\n", @commands);
# my $xspec = `which xspec`;
# chomp $xspec;
# if(!open SCR, ">$script"){
# $log->error(1, "failed to open file $script to run xspec during XrtGrbLc, exit 1");
# return;
# }
# print SCR "$str\n";
# close SCR;
# $log->entry("Running $xspec\n");
# my $xselstr = `export HEADAS=$headas; source $headas/headas-init.sh; $xspec \-$script`;
# $log->text($xselstr);
$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;
}
unlink @localclean;
chdir $curdir;
}
1