package Subs::XrtGrbLc; ############################################################################## # # DESCRIPTION: Produce XRT GRB lightcurve gif file # # HISTORY: # HISTORY: $Log: XrtGrbLc.pm,v $ # 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 ############################################################################# # 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 # 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 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 # # 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 # 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 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' ); # 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