Subs::Spectra (version 0.0)

package Subs::Spectra;
##############################################################################
#
# DESCRIPTION: This subroutine runs extracts spectra from event data
#
# HISTORY: 
# HISTORY: $Log: Spectra.pm,v $
# HISTORY: Revision 1.4  2004/05/06 20:02:34  dah
# HISTORY: Add version number back into the header comments.
# HISTORY:
# HISTORY: Revision 1.3  2004/04/16 20:21:18  dah
# HISTORY: Begin using embedded history records
# HISTORY:
#
# VERSION: 0.0
#
#
##############################################################################


use Subs::Sub;
use Util::Xanadu;

@ISA = ("Subs::Sub");
use strict;

sub new {
    my $proto=shift;
    my $self=$proto->SUPER::new();

    $self->{DESCRIPTION}="Extracting and plotting spectra";

    return $self;
}

##################
# METHODS:
##################

sub body {
    my $self=shift;

    my $log     =$self->log();
    my $filename=$self->filename();
    my $procpar =$self->procpar();
    my $jobpar  =$self->jobpar();

    $log->entry("Extracting spectra from event data");


    ###########################
    # fire up the extractor
    ###########################
    my $extractor = Util::Extractor->new();
    $extractor->verbose(2);

    ########################################
    # loop over instruments
    ########################################
    my $inst;
    for $inst ("xrt", "uvot") {

        $log->entry("Extracting $inst spectra");

        ################################################
        # determine the modes we should look at 
        ################################################
        my @modes=();
        if($inst eq "xrt") { @modes=('pc*', 'wt*', 'lr*') }

        ########################
        # loop over modes 
        ########################
        my $mode;
        foreach $mode (@modes) {

            $log->entry("Mode $mode");
            $extractor->instrument($inst, $mode);

            ######################################
            # get the input and output file names 
            ######################################
            my $spec = $filename->get("spectrum", $inst, $mode, 001);
            unlink $spec;

            my @evts = $filename->get("event", $inst, $mode, "*");
            unless(@evts) {
                $log->entry("No event files");
                next;
            }

            my $list = Util::EventFileList->new(@evts);

            $log->entry("Extracting $spec from @evts");

            ##############################
            # run the extractor
            ##############################
            $extractor->infiles(@evts)
                      ->outfile($spec, "spectrum");

            ################################################
            # this is a cludge since XRT bright mode files
            # currently lack TLMAX for PHA/PI 2002-10-22
            # This code should be removed when the problem 
            # is fixed.
            ################################################
            my $phamax = $extractor->{PARAMS}->{phamax};
            my $fits = Util::FITSfile->new($extractor->{SAMPLE_INPUT}, 
                                           "EVENTS");

            unless(defined $fits->keyword($phamax)) {
                $log->error(1, "Event file does not have $phamax keyword");
                $log->entry("Setting $phamax to 4096");

                foreach (@evts) {
                    Util::FITSfile->new($_, "EVENTS")
                                  ->keyword($phamax, 4096);
                }
            }


            $extractor->run();

        } # end of loop over modes
    } # end of loop over instruments

    ##############################
    # plot the lightcurves
    ##############################
    $self->plot_spectra();


} # end of body method

############################################################################
# plot all the lightcurves
############################################################################
sub plot_spectra {
    my $self=shift;

    my $log     =$self->log();
    my $filename=$self->filename();
    my $procpar =$self->procpar();
    my $jobpar  =$self->jobpar();

    $log->entry("Plotting all spectra");

    #################################
    # set up xspec
    #################################
    my $xspec     = Util::Xanadu->new("xspec");

    ###############################
    # loop over all spectra
    ###############################
    my $spectrum;
    foreach $spectrum ($filename->any("spectrum") ) {

        #######################################
        # get the names of corresponding files
        #######################################
        my $specplot=$filename->corresponding("spectrum","specplot",$spectrum);
    #    my $rmf=$filename->corresponding("spectrum","rmf",$spectrum);
    #    my $arf=$filename->corresponding("spectrum","arf",$spectrum);

    #    unless( -s $rmf && -s $arf ) {
    #        #############################
    #        # some files don't exist
    #        #############################
    #        $log->entry("No $arf or $rmf, so skipping $specplot");
    #        next;
    #    }


        ###############################
        # generate the plotting script
        ###############################
        $xspec->script("data $spectrum",
                       "setplot channel",
                       "ignore bad",
                       "setplot com log y on",
                       "setplot com rescale y",
                       "setplot com hard /ps",
                       "plot",
                       "exit" );

        ###########################
        # plot the spectrum
        ###########################
        $log->entry("Plotting $specplot from $spectrum");
        $xspec->run();

        ######################################
        # give the plot file the correct name
        ######################################
        rename "pgplot.ps", $specplot;

    } # end of loop over spectra

} # end of plot_spectra method