package Subs::BATSpectra; ############################################################################## # # DESCRIPTION: Accumulate BAT spectra and follow on products # # Port of Craig Markwardt's BAT::spectra, BAT::lc and setup to SDC. # # # HISTORY: $Log: BATSpectra.pm,v $ # HISTORY: Revision 1.34 2007/02/05 15:27:58 apsop # HISTORY: Trap errors from batmasktaglc that warn of no GTIs. # HISTORY: # HISTORY: Revision 1.33 2007/01/31 21:18:18 apsop # HISTORY: Do not make products from slew-only event data. # HISTORY: # HISTORY: Revision 1.32 2006/10/26 21:27:31 apsop # HISTORY: Apply GTI_TOT when extracting spectra. Trap no good time warning. Update spectra keywords. # HISTORY: # HISTORY: Revision 1.31 2006/09/28 18:31:59 apsop # HISTORY: Fix bug in time intervals used for producing bat spectra, which were completely wrong. # HISTORY: # HISTORY: Revision 1.30 2006/07/25 19:49:34 apsop # HISTORY: Add in calls to batupdatephakw and batphasyserr. # HISTORY: # HISTORY: Revision 1.29 2006/06/28 17:12:16 apsop # HISTORY: Use new target pointing GTI for filtering mask weighted light curves. # HISTORY: # HISTORY: Revision 1.28 2006/06/15 22:21:36 apsop # HISTORY: Remove "unsettled" time mask tagged lightcurves, and remove empty lightcurves. # HISTORY: # HISTORY: Revision 1.27 2006/05/06 21:16:37 apsop # HISTORY: Set enable_map to none if one is not found. # HISTORY: # HISTORY: Revision 1.26 2005/12/19 16:13:18 apsop # HISTORY: Add date keywords to primary header of mask tagged lightcurves. # HISTORY: # HISTORY: Revision 1.25 2005/11/08 17:28:44 apsop # HISTORY: Change calls for qualcal filenames. Comment out temporary fix to enable map extension name. # HISTORY: # HISTORY: Revision 1.24 2005/08/24 12:41:45 apsop # HISTORY: Bugzilla SDC 274 [from Craig Markwardt]. # HISTORY: Summary: BAT GRB processing - make RATE spectra instead of COUNTS # HISTORY: "The SDC pipeline has been defaulting to making COUNTS spectra # HISTORY: for BAT GRB processing. # HISTORY: "It turns out that the standard downstream processing tools like # HISTORY: GRPPHA and CMPPHA are not very good at handling type II counts # HISTORY: spectra. As a workaround, I recommend that the SDC switch # HISTORY: the BAT spectra to be RATE spectra." # HISTORY: # HISTORY: Revision 1.23 2005/07/05 14:56:17 apsop # HISTORY: Caught SDC up to date with BAT pipeline with exception of BAT # HISTORY: position refinement. # HISTORY: # HISTORY: Revision 1.22 2005/06/01 16:02:19 apsop # HISTORY: Change method for fetching quad rate files, as rate types have change. # HISTORY: # HISTORY: Revision 1.21 2005/05/16 14:15:33 apsop # HISTORY: Handle missing input files more gracefully when creating images. Pass # HISTORY: CALDB for ebounds parameter when creating mask tagged light curves. # HISTORY: # HISTORY: Revision 1.20 2005/05/13 19:38:02 apsop # HISTORY: Implemented product agreements between HEASARC and BAT team. # HISTORY: # HISTORY: Revision 1.19 2005/05/04 21:29:28 apsop # HISTORY: Create mask tagged light curves. # HISTORY: # HISTORY: Revision 1.18 2005/04/08 01:15:14 apsop # HISTORY: Remove unweighted 4ms lc # HISTORY: # HISTORY: Revision 1.17 2005/04/06 15:41:53 apsop # HISTORY: Change to using CALDB for cal parameters. # HISTORY: # HISTORY: Revision 1.16 2004/12/31 01:36:26 apsop # HISTORY: Comment out production of spectra which are not yet official. # HISTORY: # HISTORY: Revision 1.15 2004/12/06 01:00:12 apsop # HISTORY: Test for presence of event files, change some error() calls than entry() calls. # HISTORY: # HISTORY: Revision 1.14 2004/11/29 22:14:35 apsop # HISTORY: Corrected efile+escale parameters to batdrmgen. Only process existing # HISTORY: files. # HISTORY: # HISTORY: Revision 1.13 2004/11/16 14:30:25 apsop # HISTORY: Avoid using undefined value. # HISTORY: # HISTORY: Revision 1.12 2004/11/09 22:22:29 apsop # HISTORY: Corrected name of peak GTI. # HISTORY: # HISTORY: Revision 1.11 2004/11/09 22:00:42 apsop # HISTORY: Collect BAT GTIs in a single file. # HISTORY: # HISTORY: Revision 1.10 2004/11/08 18:41:54 apsop # HISTORY: Implemented Craig's Markwardt's BAT::spectra and BAT::lc modules. # HISTORY: # # VERSION: $Revision: 1.34 $ # # ############################################################################## use Subs::Sub; @ISA = ("Subs::Sub"); use strict; use Util::HEAdas; use Util::BATCave; sub new { my $proto=shift; my $self=$proto->SUPER::new(); $self->{DESCRIPTION}= 'Build BAT spectra, response matrices, light curves'; return $self; } ################## # METHODS: ################## sub body { my $self = shift; my $log = $self->log; my $filename = $self->filename; my $procpar = $self->procpar; my $jobpar = $self->jobpar; ########################################################## # Bail if this is slew only data not from a burst ########################################################## my @eventFiles = $filename->getExisting('unfiltered', 'bat', 'evsh{ps,sl,as}', '*'); my $seq = $jobpar->read('sequence'); return if( $seq%1000 != 0 && $seq%1000 < 899 && @eventFiles == 1 && $eventFiles[0] =~ /evshsl/ ); $self->make_spectra; $self->make_response_matrices; $self->make_light_curves; $self->make_mask_tagged_light_curves; } sub make_spectra { my $self = shift; my $log = $self->log; my $filename = $self->filename; my $procpar = $self->procpar; my $jobpar = $self->jobpar; $log->entry('Accumulating spectra'); my @events; my %info = ( ebins => 'CALDB', gti => $filename->get('burstgti', 'bat', '', 0) . '[GTI_TOT]', files => \@events, qmap => $filename->get('qualcal', 'bat', ''), auxfile => $filename->get('eventaux', 'b', '', 0) ); ########################################################### # only the preslew, slew and postslew spectra are HEASARC products ########################################################### @events = $filename->get('unfiltered', 'bat', 'evshps', '*'); if( @events ){ $info{outfile} = $filename->get('spectrum', 'bat', 'evps', 0); $self->makeSpectrum(\%info) } @events = $filename->get('unfiltered', 'bat', 'evshas', '*'); if( @events ){ $info{outfile} = $filename->get('spectrum', 'bat', 'evas', 0); $self->makeSpectrum(\%info) } @events = $filename->get('unfiltered', 'bat', 'evshsl', '*'); if( @events ){ $info{outfile} = $filename->get('spectrum', 'bat', 'evsl', 0); $self->makeSpectrum(\%info) } } sub makeSpectrum { my ($self, $info) = @_; my $log = $self->log; my $inlist = join(',', @{ $info->{files} }); my $qmap = $info->{qmap} || 'NONE'; my $gti = $info->{gti} || 'NONE'; unlink($info->{outfile}); ############################################## # set up batbinevt ############################################# my $binner = Util::HEAdas->new('batbinevt') ->params({ infile => $inlist, outfile => $info->{outfile}, outtype => 'PHA', timedel => 0, timebinalg => 'uniform', energybins => $info->{ebins}, detmask => $qmap, ecol => 'ENERGY', weighted => 'YES', outunits => 'RATE', gtifile => $info->{gti}, clobber => 'YES', tstart => 'INDEF', tstop => 'INDEF', snrthresh => 6.0, buffersize => 16384, chatter => 5, history => 'YES', }); $binner->seriousness(1); $binner->run(); if( $binner->had_error() == 2 ){ $log->error(2, "batbinevt level 2 error"); }elsif( $binner->had_error() == 1 ){ ################################################################## # Note: this case is not fatal, *IF* there were no overlapping # good times between the input file and the good time interval # file. In that case batbinevt does not create an output file. ################################################################# my $errout = $binner->stderr(); foreach ( $errout =~ /^.*$/gm ){ if( ! /^=+$/ && ! /WARNING: no overlapping good time intervals were found/ ){ $log->error(2, "batbinevt level 2 error"); last; } } } #################### # check for errors #################### if ($binner->had_error || ! -s $info->{outfile}) { $log->entry("Did not sucessfully make a spectrum for $inlist"); return; } my $specfits = Util::FITSfile->new($info->{outfile}, 1); my ($tstart, $tstop) = ( $specfits->keyword('TSTART'), $specfits->keyword('TSTOP') ); my $exp = $specfits->keyword('EXPOSURE'); my $start = Util::Date->new($tstart); my $stop = Util::Date->new($tstop); my $nhdus = $specfits->nhdus(); for(my $ext=0; $ext<$nhdus; $ext++){ next if $ext == 1; $specfits->ext($ext); $specfits->begin_many_keywords(); $specfits->keyword('TSTART', $tstart ); $specfits->keyword('TSTOP', $tstop ); $specfits->keyword('DATE-OBS', $start->date() .'T'.$start->time() ); $specfits->keyword('DATE-END', $stop->date() .'T'.$stop->time() ); $specfits->keyword('EXPOSURE', $exp ); $specfits->end_many_keywords(); } if( -f $info->{auxfile} ){ Util::HEAdas->new('batupdatephakw') ->params({infile => $info->{outfile}, auxfile => $info->{auxfile} }) ->run(); } Util::HEAdas->new('batphasyserr') ->params({infile => $info->{outfile}, syserrfile => 'CALDB' }) ->run(); } sub make_response_matrices { my $self = shift; my $log = $self->log; my $filename = $self->filename; my $procpar = $self->procpar; my $jobpar = $self->jobpar; $log->entry('Computing response matrices'); my $hkfile = $filename->get('bdaphk', 'bat') || 'NONE'; my %common = ( calfile => 'CALDB', depthfile => 'CALDB', escale => 'FILE', efile => 'CALDB', hkfile => $hkfile, chatter => 3, clobber => 'yes', ); my $batdrmgen = Util::HEAdas->new('batdrmgen') ->params(\%common); # preslew my $spectrum = $filename->get('spectrum', 'bat', 'evps'); if (-f $spectrum) { $batdrmgen->params({ infile => $spectrum, outfile => $filename->get('response', 'bat', 'evps', 0), }) ->run; } # postslew $spectrum = $filename->get('spectrum', 'bat', 'evas'); if (-f $spectrum) { $batdrmgen->params({ infile => $spectrum, outfile => $filename->get('response', 'bat', 'evas', 0), }) ->run; } } sub make_light_curves { my $self = shift; my $log = $self->log; my $filename = $self->filename; my $procpar = $self->procpar; my $jobpar = $self->jobpar; $log->entry('Accumulating light curves'); my @events = $filename->getExisting('unfiltered', 'bat', 'evsh{ps,sl,as}'); if (not @events) { $log->entry('no BAT preslew,slew,postslew events, not creating light curves'); return; } my @config = ( { desc => 'weighted 1s light curve', outfile => $filename->get('lightcurve', 'bat', 'ev1s', 0), ebins => Util::BATCave::chan4, gti => 'NONE', tbinsize => 1.0, weighted => 'YES', }, { desc => 'weighted 64ms 1-channel light curve', outfile => $filename->get('lightcurve', 'bat', 'evms', 0), ebins => Util::BATCave::chan4, gti => 'NONE', tbinsize => 0.064, weighted => 'YES', }, ); my $qmap = $filename->get('qualcal', 'bat', '') || 'NONE'; foreach my $c (@config) { $c->{files} = \@events; $c->{qmap} = $qmap; $self->lcMake($c); } } # BAT::lc->make sub lcMake { my ($self, $info) = @_; my $log = $self->log; my $filename = $self->filename; my $procpar = $self->procpar; my $jobpar = $self->jobpar; my %par = ( infile => join(',', @{ $info->{files} }), outfile => $info->{outfile}, outtype => 'LC', energybins => $info->{ebins}, detmask => $info->{qmap}, weighted => $info->{weighted} || 'NO', outunits => 'RATE', ecol => 'ENERGY', gtifile => $info->{gti} || 'NONE', clobber => 'yes', ); if ($info->{tbinsize} =~ /^g/i) { $par{timedel} = 0; $par{timebinalg} = 'gti'; } else { $par{timedel} = $info->{tbinsize}; $par{timebinalg} = 'uniform'; } my $tool = Util::HEAdas->new('batbinevt') ->params(\%par) ->run; } sub make_mask_tagged_light_curves { my $self = shift; my $log = $self->log; my $filename = $self->filename; my $procpar = $self->procpar; my $jobpar = $self->jobpar; $log->entry('Making mask tagged light curves'); my $enable_map = undef; my $tool = Util::HEAdas->new('batmasktaglc') ->params({ ebounds => 'CALDB', }); $tool->seriousness(1); my $pointing = $filename->get('gti', 's', 'tp', 0); my $tempwt = 'maskwt.tmp'; foreach my $infile ($filename->get('rawmtlc', 'bat')) { ########################################################### # Filter out times when the s/c wasn't settled, and remove # any resulting empty files ########################################################### if( -f $pointing ){ my $copy = Util::HEAdas->new('ftcopy') ->params({infile => $infile . '[1][gtifilter(\"'. $pointing .'\")]', outfile => $tempwt}); $copy->run(); my $newwt = Util::FITSfile->new($tempwt); if($newwt->nrows() == 0){ $log->entry("Deleting $infile - empty after filtering."); unlink $infile, $tempwt; next; } my $stats=Util::HEAdas->new("ftstat") ->params({infile => $tempwt .'[1][col TIME]', centroid => 'no'}) ->verbose(0); $stats->run(); my $parfile = $stats->parfile(); my ($pstart, $pstop) = ( $parfile->read('min'), $parfile->read('max') ); my $start = Util::Date->new( $pstart ); my $stop = Util::Date->new( $pstop ); foreach my $ext (0,1){ $newwt->ext($ext); $newwt->begin_many_keywords(); $newwt->keyword('TSTART', $pstart ); $newwt->keyword('TSTOP', $pstop ); $newwt->keyword('DATE-OBS', $start->date() .'T'.$start->time() ); $newwt->keyword('DATE-END', $stop->date() .'T'.$stop->time() ); $newwt->end_many_keywords(); } unlink $infile; rename $tempwt, $infile; } my $start = 0; my $status = Util::SimpleFITS->open('<' . $infile) ->readkey(TSTART => $start) ->move('MASK_TAG_WEIGHT') ->close ->status; if ($status) { $log->entry("$infile does not have MASK_TAG_WEIGHTs"); next; } $enable_map = $filename->fetch_from_repository('bdetflag', 'bat', '', $start); if (not $enable_map) { $log->error(1, "unable to get BAT enable map for time $start"); $enable_map = 'none'; } my $outfile = $infile; $outfile =~ s/_rw//; unlink($outfile) or $log->error(1, "unable to remove $outfile [$!]") if -e $outfile; my $quadfile = (grep /brtqd\.lc/, $filename->get('batrates', 'b', ''))[0]; if (not $quadfile) { $log->error(1, "missing quadrate file for $infile"); next; } $tool->params({ infile => $infile, quadfile => $quadfile, maskwt => $infile, # BAT2FITS appends the mask weight table to the raw LC outfile => $outfile, detmask => $enable_map, }) ->run(); if( $tool->had_error() == 2 ){ $log->error(2, "batmasktaglc level 2 error"); }elsif( $tool->had_error() == 1 ){ my $errout = $tool->stderr(); foreach ( $errout =~ /^.*$/gm ){ unless( /WARNING: There were no GTIs. Is that right?/ ){ $log->error(2, "batmasktaglc level 2 error"); last; } } } ####################################### # Copy DATE keywords to primary header ####################################### my $mtfits = Util::FITSfile->new($outfile); my $obs = $mtfits->keyword('DATE-OBS'); my $end = $mtfits->keyword('DATE-END'); $mtfits->ext(0); $mtfits->keyword('DATE-OBS', $obs); $mtfits->keyword('DATE-END', $end); } } 1;