astro: Add sfr_chandra.pl (version: 2015/03/04)

This script is used to do Chandra automatic spectrum extraction and
fitting.  It was written by Junhua GU and has been significantly updated
by me.

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: #