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