Util::Xspec (version 1.2)

package Util::Xspec;
##############################################################################
#
# DESCRIPTION: This Xanadu subclass adds functionality specific to xspec
# DESCRIPTION: specifically a method to extract the results of model fits.
#
# HISTORY
# HISTORY: $Log: Xspec.pm,v $
# HISTORY: Revision 1.2  2006/08/01 20:35:34  apsop
# HISTORY: Add in CVS history indicator.
# HISTORY:
# HISTORY: 1.0 -> 1.1 2002-04-01
# HISTORY: Now allows an Xserver display to be specified when initializing
# HISTORY:
# HISTORY: 1.1 -> 1.2 2003-05-16
# HISTORY: Added more methods for parsing chi-squared and flux from stdout.
# 
# VERSION: 1.2
#
##############################################################################

use Util::Xanadu;
@ISA=("Util::Xanadu");
use strict;

sub new { 
    my $self=shift;

    return $self->SUPER::new("xspec");   
}


##############################################################################
# parse the results of a model fit.
# returns a list of hash references. There is one hash per model parameter.
# The keys to the hash are as follows:
# FIT - the number of the model fit - more than one "fit" command can be given
#       in a single xspec session.
# COMPONENT  - the model component number as assigned by xspec
# MODEL      - the name of the xspec model for which this is a component.
# PARAM      - the name of the parameter within the model
# VALUE      - the best fit value
# ERROR      - an estimate for the error in the best fit value. Note that
#              these values are estimated from the second derivative of 
#              chi-squared at the minimum and are not as reliable as errors
#              derived from an exploration of chi-squared space around the 
#              minimum.
# 
##############################################################################
sub fit_results {
    my $self = shift;

    ############################
    # get the stdout text
    ############################
    my $text = $self->stdout();

    ###############################################
    # split by commands and their following text
    ###############################################
    my @commands = split /^[\s!]*XSPEC>\s*/m, $text;

    #############################
    # loop over all fit commands
    #############################
    my @results=();
    my $fit_number=0;
    my $fit;
    foreach $fit (grep /^fit/, @commands) {

        #############################################################
        # pull out the fit results these are offset by dashed lines
        #############################################################
        my ($result) = $fit =~ /^\s*-+\s*^\s*-+\s*^(.*)^\s*-+\s*^\s*-+\s*$/sm; 

        #########################################################
        # split the results into lines and select the numbered
        # lines since these are the parameters
        #########################################################
        my @lines = split /\n/, $result;
        @lines = grep /^\s*[1-9]/, @lines;

        ###################################
        # loop over lines parsing each one 
        ###################################
        foreach (@lines) {

            my @fields = split;

            my $record = {};
            $record->{FIT}      = $fit_number;
            $record->{COMPONENT} = $fields[2];
            $record->{MODEL}    = $fields[3];
            $record->{PARAM}    = $fields[4];
            $record->{VALUE}    = $fields[@fields-3];
            $record->{ERROR}    = $fields[@fields-1];

            push @results, ($record);

        } # end of loop over parameters

        $fit_number++;       
    } # end of loop over fits

    return @results;

} # end of fit_results method

##############################################################################
#
##############################################################################
sub chi_squared {
    my $self = shift;

    ############################
    # get the stdout text
    ############################
    my $text = $self->stdout();

    ###############################################
    # split by commands and their following text
    ###############################################
    my @commands = split /^[\s!]*XSPEC>\s*/m, $text;

    #############################
    # loop over all fit commands
    #############################
    my @results=();
    my $fit_number=0;
    my $fit;
    foreach $fit (grep /^fit/, @commands) {

        my $record={};
        $record->{FIT} = $fit_number++;

        my ($chi2)    = $fit =~                 /Chi-Squared =\s+([^\s]+)/;
        my ($reduced) = $fit =~         /Reduced chi-squared =\s+([^\s]+)/;
        my ($prob)    = $fit =~ /Null hypothesis probability =\s+([^\s]+)/;
        my ($dof)     = $fit =~ /(\d+) degrees of freedom/;

        $record->{CHI2}        = $chi2;
        $record->{REDUCED_CHI2}= $reduced;
        $record->{PROBABILITY} = $prob;
        $record->{FREEDOM}     = $dof;


        #print "chi2=$chi2 reduced=$reduced dof=$dof prob=$prob\n";
        
        push @results, ($record);

    } # end of loop over fit commands
  
    return @results;

} # end of chi_squared method

##############################################################################
#
##############################################################################
sub flux {
    my $self = shift;

    ############################
    # get the stdout text
    ############################
    my $text = $self->stdout();

    ###############################################
    # split by commands and their following text
    ###############################################
    my @commands = split /^[\s!]*XSPEC>\s*/m, $text;

    #############################
    # loop over all fit commands
    #############################
    my @results=();
    my $fit_number=0;
    foreach my $fit (grep /^flux/, @commands) {

        my $record={};

        my ($photons) = $fit =~ /([^\s]+)\s+photons/;
        my ($ergs) = $fit =~ /([^\s]+)\s+ergs/;

        my @errors = $fit =~ 
   /Error range\s+([^\s]+)\s+-\s+([^\s]+)\s+\(\s*([^\s]+)\s+-\s+([^\s]+)\s*\)/;

        $record->{FIT}     = $fit_number++;
        $record->{PHOTONS} = $photons;
        $record->{ERGS}    = $ergs;

        $record->{PHOTONS_LO} = shift @errors;
        $record->{PHOTONS_HI} = shift @errors;
        $record->{ERGS_LO}    = shift @errors;
        $record->{ERGS_HI}    = shift @errors;


        push @results, ($record);

    }

    return @results;

} # end of flux command


sub DESTROY {
    my $self = shift;
    # check for an overridden destructor...
    $self->SUPER::DESTROY if $self->can("SUPER::DESTROY");
}


1;