diff options
Diffstat (limited to 'astro/spectrum')
| -rwxr-xr-x | astro/spectrum/sfr_chandra.pl | 686 |
1 files changed, 686 insertions, 0 deletions
diff --git a/astro/spectrum/sfr_chandra.pl b/astro/spectrum/sfr_chandra.pl new file mode 100755 index 0000000..a6f982c --- /dev/null +++ b/astro/spectrum/sfr_chandra.pl @@ -0,0 +1,686 @@ +#!/usr/bin/env perl +# +use strict; +use warnings; +use File::Basename; +# +# Performs automatically spectral fitting with XSPEC 12. +# For Chandra data only. +# +# Reference: +# [1] Running Multiple Copies of a Tool - CIAO +# http://cxc.harvard.edu/ciao/ahelp/parameter.html#Running_Multiple_Copies_of_a_Tool +# +# Author: Junhua GU <tompkins@sjtu.edu.cn> +# 2007 +# +# Updated by: Weitian LI <liweitianux@gmail.com> +my $version = "2015/03/04"; +# +# ChangeLogs: +# 2015/03/04: Weitian LI +# Totally rewrite of this script: +# * add use strict, warnings, many fixes to the syntax +# * add sub 'usage' +# * rewrite sub 'update_backscale' -> 'rescale' +# * update sub 'extract_spec' -> 'extract_spectrum' +# * remove sub 'extract_back' +# * replace 'grppha' with 'dmgroup' +# * add sub 'determine_bkg' to determine the type of background +# * fix the problem of when background set to 'none' +# * copy CIAO tools' parameter files to local, avoiding conflicts when +# running multiple instances of a CIAO tool +# * allow comments in fit config file (starting with '//') +# * update sub 'parse_result' to accept the parameter list and hash, +# for better/simple parse process +# * update example fit config file +# * update generation of xspec fit script +# * update config file syntax: merged '#PARA' and '#ERROR' +# * use backtics to get the output of command within perl +# * !! INCOMPATIBLE !! improve/change output format (columns adjustments) +# 2015/02/05: Weitian LI +# * fix the wrong use of chomp. +# 2015/01/27: Weitian LI +# * check 'cnt_spec' value before adjust the backscale +# 2015/01/21: +# * add 'query yes' to temp_fit_script +# + +## usage {{{ +sub usage { + my $prog = basename($0); + my $text = qq(SFR - Spectral Fitting Robot +Usage: + $prog <region_list> <evt> <bkgd> <output> <rmf> <arf> <fit_config> + +Version: $version +); + print $text; +} +## usage }}} + +# Example fit config {{{ +#---------------------------------------------------------------- +# abund grsa +# weight gehrels +# +# model wabs*apec & <nh> & 1.0 & 0.5 & <redshift> & 1.0 & +# +# freeze 1 +# thaw 3 +# +# ignore **:0.0-0.7,7.0-** +# +# fit +# fit +# fit +# +# // Do not adjust BACKSCAL of background spectrum +# #NORESCALE +# +# #PARA kT 2 +# +# // Also calculate the errors of kT with XSPEC error command +# #PARA kT 2 ERROR 1.0 +# #PARA Abund 3 +# #STAT +# #FLUX 0.7 7.0 +# +# // Specify the parameters for dmgroup +# #DMGROUP NUM_CTS 20 +# #DMGROUP BIN 1:128:2,129:256:4,257:512:8,513:1024:16 +#---------------------------------------------------------------- +# example config }}} + +## variables {{{ +# Background types +use constant { + NO_BKG => 0, # do not use a background + BLANKSKY => 1, # blanksky fits file + STOW_BKG => 2, # stowed background + LBKG_REG => 3, # region file of a local background + BKG_SPEC => 4, # background spectrum file +}; + +# Threshold value for the 'reduced chi-squared' +my $rechisq_threshold = 1.4; +# Number of bad fitting (If fitted rechisq > $rechisq_threshold) +my $bad_fit_num = 0; + +# Whether to rescale the background (default yes) +my $rescale = 1; + +# Energy range selected to adjust the BACKSCAL of the background spectra. +# PI = [ energy(eV) / 14.6 eV + 1 ] +# Chandra particle background energy range: +# 9.5-12.0 keV (channel: 651-822) +my $pb_chn_lo = 651; +my $pb_chn_hi = 822; + +# Source spectrum particle background counts ($cnt_spec, +# within channel back_lo -> back_hi) may be too small, even ZERO, +# thus cause error of zero division. +# If '$cnt_spec' too small (say 20), then skip adjusting the BACKSCAL. +my $pb_cnt_threshold = 20; + +# 'dmgroup' parameters +my $grouptype = "NUM_CTS"; +#my $grouptype="BIN" +my $grouptypeval = "20"; +my $binspec = "1:128:2,129:256:4,257:512:8,513:1024:16"; +## variables }}} + +# NOTE: @ARGV does NOT include $0, the name of the program. +#my $argc = @ARGV; +my $argc = $#ARGV + 1; +#print "argc: $argc\n"; +if ($argc != 7) { + usage(); + exit 1; +} + +## Setup local PFILES environment {{{ +# Allow running multiple copies of a CIAO tool. +# Reference: http://cxc.harvard.edu/ciao/ahelp/parameter.html#Running_Multiple_Copies_of_a_Tool +# +# Make a local copy of necessary parameter files. +print "Setup local minimal environment for CIAO to avoid pfile conflicts.\n"; +print "Copy necessary parameter files ...\n"; +my @ciao_tools = ("dmstat", "dmkeypar", "dmhedit", "dmextract", "dmgroup"); +foreach (@ciao_tools) { + my $tool = $_; + my $tool_par = `paccess $tool`; + chomp($tool_par); + #print "${tool}.par: $tool_par.\n"; + system("punlearn $tool && cp -Lv $tool_par ."); +} + +# Set environment variable PFILES +$ENV{PFILES} = "./:$ENV{PFILES}"; +#system("echo PFILES: \$PFILES"); +## PFILES }}} + +## process arguments {{{ +# File contains the list of fitting regions +my $region_list_fn = $ARGV[0]; +# Event file (clean) +my $evt2_fn = $ARGV[1]; +# Background file (blanksky.fits / bkg.pi) +my $bkg_fn = $ARGV[2]; +# Output file to save fitting results +my $output_fn = $ARGV[3]; +my $rmf_fn = $ARGV[4]; +my $arf_fn = $ARGV[5]; +# Fit config file for this tool +my $config_fn = $ARGV[6]; + +print "====================================================\n"; +print "region_list: $region_list_fn\n"; +print "evt2: $evt2_fn\n"; +print "bkg: $bkg_fn \n"; +print "output: $output_fn\n"; +print "rmf: $rmf_fn\n"; +print "arf: $arf_fn\n"; +print "config: $config_fn\n"; + +# Parse evt, back, rmf, arf ... {{{ +my @evt_list = split(/[,\ ]+/, $evt2_fn); +my $evt_num = @evt_list; +my @bkg_list = split(/[,\ ]+/, $bkg_fn); +my $bkg_num = @bkg_list; +my @rmf_list = split(/[,\ ]+/, $rmf_fn); +my $rmf_num = @rmf_list; +my @arf_list = split(/[,\ ]+/, $arf_fn); +my $arf_num = @arf_list; + +print "====================================================\n"; +print "evt_num: $evt_num, bkg_num: $bkg_num, rmf_num: $rmf_num, arf_num: $arf_num\n"; +if (! (($evt_num == $bkg_num) && ($bkg_num == $rmf_num) && ($rmf_num == $arf_num))) { + print "ERROR: number of evt, bkg, rmf or arf NOT match!\n"; + exit 12; +} + +foreach my $i (0 .. $#evt_list) { + printf "Data group %d: ", ($i+1); + printf "%s, %s, %s, %s\n", $evt_list[$i], $bkg_list[$i], $rmf_list[$i], $arf_list[$i]; +} +print "====================================================\n"; +# parse evt, bkg ... }}} + +# Backup the output file if it exists +if ( -e $output_fn ) { + print "WARNING: Output file '$output_fn' already exists!\n"; + system("mv -fv $output_fn ${output_fn}_bak"); +} + +# Determine the type of each background file +my @bkg_type = (); +foreach my $bkgd (@bkg_list) { + my $bkgd_type = NO_BKG; + if ($bkgd =~ /^(no|null|none)$/i) { + print "WARNING: does not use background!\n"; + $bkgd_type = NO_BKG; + } elsif (! -e $bkgd) { + print "ERROR: bkg file '$bkgd' does not exists!\n"; + exit 11; + } else { + $bkgd_type = determine_bkg($bkgd); + } + push(@bkg_type, $bkgd_type); +} +## arguments }}} + +## Fitting related variables {{{ +my $fit_reg_fn = "_fit.reg"; +my $fit_script_fn = "_fit.tcl"; +my $fit_result_fn = "_fit.txt"; +if (-e $fit_result_fn) { + unlink($fit_result_fn); +} + +my @fit_spec_list = (); +my @fit_grp_spec_list = (); +my @fit_bkg_spec_list = (); +foreach my $i (0 .. $#evt_list) { + my $ii = $i + 1; + push (@fit_spec_list, "_fit_${ii}.pi"); + push (@fit_grp_spec_list, "_fit_${ii}_grp.pi"); + push (@fit_bkg_spec_list, "_fit_${ii}_bkg.pi"); +} +## varaibles }}} + +## Subroutines {{{ +# Determine the type of the given background file {{{ +# Usage: determine_bkg($bkg_filename) +sub determine_bkg { + my $argc = @_; + if ($argc != 1) { + print "ERROR: sub determine_bkg() requires 1 argument!\n"; + exit 21; + } + my $bkgd = $_[0]; + my $filetype = `file -bL $bkgd`; + if ($filetype =~ /^ASCII text/) { + print "determine_bkg: '$bkgd' is ASCII text\n"; + return LBKG_REG; + } elsif ($filetype =~ /^FITS/) { + print "determine_bkg: '$bkgd' is FITS\n"; + my $hduclas1 = `punlearn dmkeypar && dmkeypar $bkgd HDUCLAS1 echo=yes`; + chomp($hduclas1); + if ($hduclas1 =~ /^EVENTS$/) { + print "determine_bkg: '$bkgd' is FITS of EVENTS\n"; + my $bkg_obj = `punlearn dmkeypar && dmkeypar $bkgd OBJECT echo=yes`; + chomp($bkg_obj); + if ($bkg_obj =~ /^BACKGROUND DATASET$/) { + print "determine_bkg: '$bkgd' is background dataset\n"; + return BLANKSKY; + } elsif ($bkg_obj =~ /^ACIS STOWED$/) { + print "determine_bkg: '$bkgd' is stowed background\n"; + return STOW_BKG; + } else { + print "ERROR: determine_bkg: unknown EVENTS type\n"; + exit 22; + } + } elsif ($hduclas1 =~ /^SPECTRUM$/) { + print "determine_bkg: '$bkgd' is FITS of SPECTRUM\n"; + return BKG_SPEC; + } else { + print "ERROR: determine_bkg: unknown FITS type\n"; + exit 23; + } + } else { + print "ERROR: determine_bkg: unknown file type: $filetype\n"; + exit 24; + } +} +# determine_bkg }}} + +# Adjust BACKSCALE of background spectrum to match PB fluxes {{{ +# Usage: bkg_rescale($src_spec, $bkg_spec) +sub bkg_rescale { + my $argc = @_; + if ($argc != 2) { + print "ERROR: sub bkg_rescale() requires 2 arguments!\n"; + exit 22; + } + my $src_spec = $_[0]; + my $bkg_spec = $_[1]; + + # PB counts of source & background spectrum + my $src_pb_cnt = `punlearn dmstat && dmstat infile="$src_spec\[channel=$pb_chn_lo:$pb_chn_hi\]\[cols COUNTS\]" >/dev/null 2>&1 && pget dmstat out_sum`; + chomp($src_pb_cnt); + # PB counts of background spectrum + my $bkg_pb_cnt = `punlearn dmstat && dmstat infile="$bkg_spec\[channel=$pb_chn_lo:$pb_chn_hi\]\[cols COUNTS\]" >/dev/null 2>&1 && pget dmstat out_sum`; + chomp($bkg_pb_cnt); + #print "src_pb_cnt: $src_pb_cnt\n"; + #print "bkg_pb_cnt: $bkg_pb_cnt\n"; + #exit 233; + + # BACKSCAL & EXPOSURE of source spectrum + my $src_backscal = `punlearn dmkeypar && dmkeypar $src_spec BACKSCAL echo=yes`; + chomp($src_backscal); + my $src_exposure = `punlearn dmkeypar && dmkeypar $src_spec EXPOSURE echo=yes`; + chomp($src_exposure); + # BACKSCAL & EXPOSURE of background spectrum + my $bkg_backscal = `punlearn dmkeypar && dmkeypar $bkg_spec BACKSCAL echo=yes`; + chomp($bkg_backscal); + my $bkg_exposure = `punlearn dmkeypar && dmkeypar $bkg_spec EXPOSURE echo=yes`; + chomp($bkg_exposure); + #print "src_backscal: $src_backscal\n"; + #print "src_exposure: $src_exposure\n"; + #print "bkg_backscal: $bkg_backscal\n"; + #print "bkg_exposure: $bkg_exposure\n"; + + if ($src_pb_cnt <= $pb_cnt_threshold) { + # Source spectrum particle background counts too small. + print "WARNING: too small counts of source spectrum within 9.5-12 keV: $src_pb_cnt!\n"; + print "Skipped BACKSCAL adjustment of background spectrum!\n"; + } else { + # Adjust BACKSCALE of background spectrum to match PB flux: + # src_pb_cnt / src_exposure / src_backscal = bkg_pb_cnt / bkg_exposure / bkg_backscal + my $bkg_backscal_new = $src_backscal * $src_exposure * $bkg_pb_cnt + / ($src_pb_cnt * $bkg_exposure); + my $cmd_str = "punlearn dmhedit && " + . "dmhedit infile=\"$bkg_spec\" filelist=none operation=add " + . "key=BACKSCAL value=$bkg_backscal_new " + . "comment='old_value: $bkg_backscal'"; + system($cmd_str); + print "Updated BACKSCAL: $bkg_backscal -> $bkg_backscal_new\n"; + } +} +# rescale }}} + +# Extract spectrum from fits {{{ +# Usage: extract_spectrum($input_fits, $output_spec, $region_file) +sub extract_spectrum { + my $argc = @_; + if ($argc != 3) { + print "ERROR: sub extract_spec() requires 3 arguments!\n"; + exit 31; + } + my $infile = $_[0]; + my $outfile = $_[1]; + my $regfile = $_[2]; + my $cmd_str = "punlearn dmextract && " + . "dmextract infile=\"$infile\[sky=region\($regfile\)\]\[bin pi\]\" " + . "outfile=\"$outfile\" wmap=\"\[bin det=8\]\" clobber=yes"; + system($cmd_str); +} +# extract spectrum }}} + +# group spectrum with 'dmgroup' {{{ +# Usage: group_spectrum($in_spec, $grp_spec, $grouptype, $grouptypeval, $binspec) +sub group_spectrum { + my $argc = @_; + if ($argc != 5) { + print "ERROR: sub group_spec() requires 5 arguments!\n"; + exit 41; + } + my $infile = $_[0]; + my $outfile = $_[1]; + my $grp_type = $_[2]; + my $grp_type_val = $_[3]; + my $bin_spec = $_[4]; + + my $cmd_str = "punlearn dmgroup && " + . "dmgroup infile=\"$infile\" outfile=\"$outfile\" " + . "grouptype=\"$grp_type\" grouptypeval=$grp_type_val " + . "binspec=\"$bin_spec\" xcolumn=CHANNEL ycolumn=COUNTS clobber=yes"; + #print $cmd_str; + system($cmd_str); +} +# group spectrum }}} + +# Parse the fitted result line {{{ +# Usage: parse_result($result_line, @para_name_list, %para_num) +sub parse_result ($\@\%) { + my $argc = @_; + if ($argc != 3) { + print "ERROR: sub parse_result() requires 3 argument!\n"; + exit 51; + } + my $line=$_[0]; + my @para_name_list = @{$_[1]}; + my %para_num= %{$_[2]}; + + my @parsed_results = (); + + foreach my $para_name (@para_name_list) { + # Number of values for this parameter + my $num = $para_num{$para_name}; + if ($line =~ /^$para_name:/) { + my @items = split(/\s+/, $line); + foreach my $i (0 .. $num) { + # items[0] is the '$para_name:' + push(@parsed_results, $items[$i]); + } + if ($para_name =~ /red.*chi.*sq/i) { + my $current_rechisq = $items[1]; + if ($current_rechisq > $rechisq_threshold) { + print "WARNING: bad quality of spectral fittings!\n"; + $bad_fit_num ++; + } + } + } + } + + return @parsed_results; +} +# parse result }}} +## subroutines }}} + +## Generate the fit script for XSPEC {{{ +if (! -e $config_fn) { + print "ERROR: config file '$config_fn' does not exist!\n"; + exit 14; +} +open(my $FIT_SCRIPT, ">", $fit_script_fn) + or die "ERROR: Could not open file '$fit_script_fn', $!"; + +# Part of load data to XSPEC +foreach my $i (0 .. $#evt_list) { + my $ii = $i + 1; + print $FIT_SCRIPT "data ${ii}:${ii} $fit_grp_spec_list[$i]\n"; + print $FIT_SCRIPT "response 1:${ii} $rmf_list[$i]\n"; + print $FIT_SCRIPT "arf 1:${ii} $arf_list[$i]\n"; + print $FIT_SCRIPT "backgrnd ${ii} $fit_bkg_spec_list[$i]\n"; +} + +# Automatically answer XSPEC query with 'yes' +print $FIT_SCRIPT "query yes\n"; +print $FIT_SCRIPT "set ff \[ open $fit_result_fn w \]\n"; + +open(my $CONFIG, "<", $config_fn) + or die "ERROR: Counld not open file '$config_fn', $!"; +my @config_lines = <$CONFIG>; +close($CONFIG); + +# Record the name of specified parameters in the config file. +my @para_name_list = (); +# Record the number of output values of above parameters. (use hash) +my %para_num; + +# Parse the fitting config file {{{ +foreach my $line (@config_lines) { + chomp($line); + if ($line =~ /^\s*$/) { + # Blank line + print "Skipped blank line.\n"; + next; + } elsif ($line =~ /^\/\//) { + # Comment line + print "Skipped comment line.\n"; + next; + } elsif ($line =~ /^\#PARA/) { + # Define fitting parameter, together with settings error calculation. + # Config: '#PARA para_name para_index [ ERROR conf_level ]' + # Output: + # 'para_name: value sigma' + # 'para_name_err: error_lower error_upper' + my @para = split(/\s+/, $line); + my $para_name = $para[1]; + my $para_idx = $para[2]; + push(@para_name_list, $para_name); + $para_num{$para_name} = 2; + print $FIT_SCRIPT "tclout para $para_idx\n"; + print $FIT_SCRIPT "scan \$xspec_tclout \"%f\" value\n"; + print $FIT_SCRIPT "tclout sigma $para_idx\n"; + print $FIT_SCRIPT "scan \$xspec_tclout \"%f\" sigma\n"; + print $FIT_SCRIPT "puts \$ff \"$para_name: \$value \$sigma\"\n"; + # Check whether to calculate the errors + if ($line =~ /ERR(OR|)\s+\d*\.\d*\s*$/) { + # Calculate errors with XSPEC 'error' command + my $conf_level = $para[4]; + my $para_err_name = "${para_name}_err"; + push(@para_name_list, $para_err_name); + $para_num{$para_err_name} = 2; + print $FIT_SCRIPT "fit\n"; + print $FIT_SCRIPT "error $conf_level $para_idx\n"; + print $FIT_SCRIPT "tclout error $para_idx\n"; + print $FIT_SCRIPT "scan \$xspec_tclout \"%f %f\" err_lower err_upper\n"; + print $FIT_SCRIPT "puts \$ff \"${para_name}_err: \$err_lower \$err_upper\"\n"; + } + } elsif ($line =~ /^\#STAT/) { + # Whether output the statistic information + # Config: '#STAT' + # Output: + # 'stat_value statistic' + # 'dof_value dof' + # 'rechisq_value Reduced_chisq' + my $stat_name = "stat"; + push(@para_name_list, $stat_name); + $para_num{$stat_name} = 1; + my $dof_name = "dof"; + push(@para_name_list, $dof_name); + $para_num{$dof_name} = 1; + # Reduced chi-squared + my $rechisq_name = "reduced_chisq"; + push(@para_name_list, $rechisq_name); + $para_num{$rechisq_name} = 1; + print $FIT_SCRIPT "tclout $stat_name\n"; + print $FIT_SCRIPT "scan \$xspec_tclout \"%f\" stat_val\n"; + print $FIT_SCRIPT "puts \$ff \"$stat_name: \$stat_val\"\n"; + print $FIT_SCRIPT "tclout $dof_name\n"; + print $FIT_SCRIPT "scan \$xspec_tclout \"%d\" dof_val\n"; + print $FIT_SCRIPT "puts \$ff \"$dof_name: \$dof_val\"\n"; + print $FIT_SCRIPT "set rechisq_cmd \{ format \"%.5f\" \[ expr \{ \$stat_val / \$dof_val \} \] \}\n"; + print $FIT_SCRIPT "set rechisq \[ eval \$rechisq_cmd \]\n"; + print $FIT_SCRIPT "puts \$ff \"$rechisq_name: \$rechisq\"\n"; + } elsif ($line =~ /^\#FLUX/) { + # Whether output flux information + # Config: '#FLUX start_keV end_keV' + my @para = split(/\s+/, $line); + my $energy_low = $para[1]; + my $energy_high = $para[2]; + my $para_name = "flux"; + push(@para_name_list, $para_name); + $para_num{$para_name} = 1; + print $FIT_SCRIPT "flux $energy_low $energy_high\n"; + print $FIT_SCRIPT "tclout $para_name\n"; + print $FIT_SCRIPT "scan \$xspec_tclout \"%f\" flux_val\n"; + print $FIT_SCRIPT "puts \$ff \"$para_name: \$flux_val\"\n"; + } elsif ($line =~ /^\#DMGROUP/) { + # dmgroup parameters + my @para = split(/\s+/, $line); + $grouptype = $para[1]; + if ($grouptype eq "NUM_CTS") { + $grouptypeval = $para[2]; + print "dmgroup: $grouptype, $grouptypeval\n"; + } elsif ($grouptype eq "BIN") { + $binspec = $para[2]; + print "dmgroup: $grouptype, $binspec\n"; + } else { + print "ERROR: unknown grouptype '$grouptype'!\n"; + print "Only 'NUM_CTS' and 'BIN' supported for dmgroup.\n"; + exit 16; + } + } elsif ($line =~ /^\#NORESCALE/) { + $rescale = 0; + } elsif ($line =~ /^\#/) { + # Other unsupported lines starting with '#' + print "WARNING: unsupported config line:\n '$line'\n"; + next; + } else { + # XSPEC commands/lines + print $FIT_SCRIPT "$line\n"; + } +} +# Parse config file }}} + +print $FIT_SCRIPT "close \$ff\n"; +print $FIT_SCRIPT "tclexit\n"; +close($FIT_SCRIPT); +#exit 233; +## Generate XSPEC script }}} + +## Perform spectral fittings {{{ +open(my $REGION_LIST, "<", $region_list_fn) + or die "ERROR: Could not open file '$region_list_fn', $!"; +my @region_lines = <$REGION_LIST>; +close($REGION_LIST); + +open(my $OUTFILE, ">", $output_fn) + or die "ERROR: Could not open file '$output_fn', $!"; + +my $reg_total = @region_lines; +my $n = 0; + +# Loop over the region list to perform spectral fitting for each region +foreach my $line (@region_lines) { + chomp($line); + $n++; + print "### $n of $reg_total ###\n"; + + # Skip the null/none region + if ($line =~ /^\#(null|none)/i) { + print $OUTFILE "$n "; + foreach my $para_name (@para_name_list) { + # Number of values for this parameter + my $num = $para_num{$para_name}; + print $OUTFILE "${para_name}: "; + foreach my $i (1 .. $num) { + print $OUTFILE "-1 "; + } + } + print $OUTFILE "\n"; + print "Skipped the null/none region.\n"; + next; + } + + # Save the current region to a temporary region file + unlink($fit_reg_fn); + open(my $FIT_REG, ">", $fit_reg_fn) + or die "ERROR: Could not open file '$fit_reg_fn', $!"; + my @cur_regs = split(/\s+/, $line); + foreach my $reg (@cur_regs) { + print $FIT_REG "$reg\n" + } + close($FIT_REG); + + # Prepare source and background spectra {{{ + print "Prepare source and background spectra \.\.\.\n"; + foreach my $i (0 .. $#evt_list) { + # Prepare source spectrum + extract_spectrum($evt_list[$i], $fit_spec_list[$i], $fit_reg_fn); + group_spectrum($fit_spec_list[$i], $fit_grp_spec_list[$i], + $grouptype, $grouptypeval, $binspec); + # Prepare background spectrum + if ($bkg_type[$i] == NO_BKG) { + # Dose not use a background + $fit_bkg_spec_list[$i] = "none"; + } elsif (($bkg_type[$i] == BLANKSKY) || ($bkg_type[$i] == STOW_BKG)) { + extract_spectrum($bkg_list[$i], $fit_bkg_spec_list[$i], $fit_reg_fn); + if ($rescale) { + bkg_rescale($fit_spec_list[$i], $fit_bkg_spec_list[$i]); + } + } elsif ($bkg_type[$i] == LBKG_REG) { + extract_spectrum($evt_list[$i], $fit_bkg_spec_list[$i], $bkg_list[$i]); + if ($rescale) { + bkg_rescale($fit_spec_list[$i], $fit_bkg_spec_list[$i]); + } + } elsif ($bkg_type[$i] == BKG_SPEC) { + # Directly use given background spectrum + system("cp -v $bkg_list[$i] $fit_bkg_spec_list[$i]"); + if ($rescale) { + bkg_rescale($fit_spec_list[$i], $fit_bkg_spec_list[$i]); + } + } else { + print "ERROR: unsupported background type!\n"; + exit 15; + } + } + # prepare spectra }}} + + # Perform the spectral fitting with XSPEC + print "Fitting \.\.\.\n"; + unlink($fit_result_fn); + system("xspec - $fit_script_fn >/dev/null") == 0 + or die "system: xspec failed: $?"; + + # Parse fitted results and output + open(my $FIT_RESULT, "<", $fit_result_fn) + or die "ERROR: Could not open file '$fit_result_fn', $!"; + my @fitted_results = <$FIT_RESULT>; + close($FIT_RESULT); + # Array to store parsed results of current fit + my @current_results = ($n); + foreach my $line (@fitted_results) { + chomp($line); + push(@current_results, parse_result($line, @para_name_list, %para_num)); + } + print join(" ", @current_results) . "\n"; + print $OUTFILE join(" ", @current_results) . "\n"; +} # end of fitting loop + +close($OUTFILE); +## spectra fittings }}} + +if ($bad_fit_num > 0) { + my $bad_fit_percent = $bad_fit_num / $reg_total; + print "\n************************************************\n"; + print "WARNING: Reduced chi-squared of $bad_fit_num ($bad_fit_percent) regions > $rechisq_threshold!!\n"; +} + +exit 0; + +# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=perl: # |