Subs::CalSources (version 0.0)


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


use Subs::Sub;
use Util::Xspec;

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

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

    $self->{DESCRIPTION}="Extracting calibration source 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 calibration source spectra from event data");

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

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

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

        ################################################
        # determine the modes we should look at 
        ################################################
        my $nsources=0;
        my @modes=();
        if($inst eq "xrt") { 
            @modes=("pc*");
            $nsources = 4;
        }

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

            $log->entry("Mode $mode");

            #####################################
            # get a list of event files
            #####################################
            my @unfs= $filename->get('unfiltered', $inst, $mode, '*');
            unless(@unfs) {
                $log->entry('No event files');
                next;
            }

            ##################################
            # set the mode for the extractor
            ##################################
            $extractor->instrument($inst, $mode);

            ################################
            # loop over calibration sources
            ################################
            my $index;
            for($index=0; $index<$nsources; $index++) {

                ############################
                # output file name
                ############################
                my $spec = $filename->get("calspec", $inst, $mode,$index);
                unlink $spec;

                my $region=$filename->fetch_cal("regcal", $inst, "", $index+1);

                $log->entry("Extracting $spec from @unfs with $region");

                ######################################
                # set the input and output file names
                ######################################
                $extractor->infiles(@unfs)
                          ->outfile($spec, "spectrum")
                          ->region($region)
                          ->params({ecol=>"PHA"});

#$log->error(1,"extracting PI calibration spectra due to problem with PHA columns");


                ################################################
                # 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 (@unfs) {
                        Util::FITSfile->new($_, "EVENTS")
                                      ->keyword($phamax, 4096);
                    }
                }


                $extractor->run();

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

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


} # end of body method

############################################################################
# plot all the spectra
############################################################################
sub fit_and_plot {
    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::Xspec->new();

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

        my ($inst, $mode, $index) = $filename->parse($spectrum, "calspec");

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

        my $rmf=$filename->fetch_cal("diagrmf", $inst);

        ###############################
        # find the peak channel
        ###############################
#        my $peak = 2800; # this was valid for Panter data
        my $peak = 2300; # this is valid for functional test

        ###############################
        # some parameters
        ###############################
        my $lo_channel = $peak-200;
        my $hi_channel = $peak+400;
        my $peak_beta  = $peak+250;

        ###############################
        # generate the plotting script
        ###############################
        $xspec->script("query yes",
                       "data $spectrum",
                       "response $rmf",

                       "ignore bad",
                       "ignore 0-$lo_channel",
                       "ignore $hi_channel-4095",

                       "model gauss/b + gauss/b",
                       "$peak", # peak channel of alpha component
                       "10.0", # width
                       " ",    # normalization

                       "$peak_beta", # peak channel of beta component
                       "10.0", # width
                       " ",    # normalization


                       "renorm",
                       "fit",

                       "setplot channel",
                       "setplot com rescale x $lo_channel $hi_channel",
                   #    "setplot com log y on",
                   #    "setplot com rescale y",

                       "setplot com hard /cps",
                       "plot",
                       "exit" );

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

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

        ####################################
        # get the fit results
        ####################################
        my @results = $xspec->fit_results;

        $log->entry("Alpha line peak = $results[0]->{VALUE} ".
                    "+/- $results[0]->{ERROR}" );

        $log->entry("Alpha line width = $results[1]->{VALUE} ".
                    "+/- $results[1]->{ERROR}" );

        $log->entry("Beta line peak = $results[3]->{VALUE} ".
                    "+/- $results[3]->{ERROR}" );

        $log->entry("Beta line width = $results[4]->{VALUE} ".
                    "+/- $results[4]->{ERROR}" );

        my $i = $index +0;
        $jobpar->set({"xrtcal_apeak".$i      => $results[0]->{VALUE},
                      "xrtcal_apeak_err".$i  => $results[0]->{ERROR},

                      "xrtcal_awidth".$i     => $results[1]->{VALUE},
                      "xrtcal_awidth_err".$i => $results[1]->{ERROR},

                      "xrtcal_bpeak".$i      => $results[3]->{VALUE},
                      "xrtcal_bpeak_err".$i  => $results[3]->{ERROR},

                      "xrtcal_bwidth".$i     => $results[4]->{VALUE},
                      "xrtcal_bwidth_err".$i => $results[4]->{ERROR} });






    } # end of loop over spectra

} # end of plot_spectra method