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.17 2014/08/15 09:33:47 apsop # HISTORY: make_xrt_grb_lc: More informative email messages to SDC staff if can't # HISTORY: findGRB in catalogs. Log a different message if it was too soon to # HISTORY: send the email. # HISTORY: # 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.17 $ # ############################################################################## 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 $seqprocnum = $jobpar->read('seqprocnum'); my $job_title = $jobpar->read('job_title'); my $obsdate = $jobpar->read('obsdate'); my $subject = "GRB not found in catalogs: $trigid"; my $msg = "Unable to find GRB info in all catalogs for:\n\n" . "object= $object\n" . "sequence= $trigid\n" . "target= $target\n" . "processing number= $seqprocnum\n" . "observation date= $obsdate\n" . "job title= $job_title\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" . "(Email generated in XrtGrbLc::make_xrt_grb_lc)\n"; $self->sendEmail($subject, $msg, $watchers); $log->error(1, "Unable to find GRB in catalogs. Email notice sent to $watchers"); } else { $log->error(1, "Unable to find GRB in catalogs. Email notice not sent, too soon"); } 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