Subs::BATSpectra (version $)


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.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.32 $
#
#
##############################################################################


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;

	$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',
				  });

  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;
	  
    #######################################
    # 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;