astro/sfr_chandra.pl: Make RMF and ARF for each region

* Use 'specextract' to extract spectrum and make RMF and ARF
* Update arguments and usage
* Many small improvements

1 files changed, 143 insertions, 89 deletions

diff --git a/astro/spectrum/sfr_chandra.pl b/astro/spectrum/sfr_chandra.pl
index aa8cd59..7dc0cc2 100755
--- a/astro/spectrum/sfr_chandra.pl
+++ b/astro/spectrum/sfr_chandra.pl
@@ -1,6 +1,6 @@
 #!/usr/bin/env perl
 #
-# Copyright (c) 2015 Weitian LI <liweitianux@gmail.com>
+# Copyright (c) 2015,2019 Weitian LI <liweitianux@gmail.com>
 # Copyright (c) 2007 Junhua GU <tompkins@sjtu.edu.cn>
 #
 # MIT License
@@ -8,20 +8,19 @@
 # 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
-#
 
 use strict;
 use warnings;
 use File::Basename;
 use File::Copy;
 
-my $version = "2019/10/16";
+my $version = "2019/10/21";
 #
 # ChangeLogs:
 #
+# 2019/10/21: Weitian LI
+#   * generate RMF and ARF for each region
+#
 # 2019/10/16: Weitian LI
 #   * support region list file with each line representing a region file
 #     instead of the region itself (Sanders's contour binning
@@ -59,22 +58,6 @@ my $version = "2019/10/16";
 #   * 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>
-
-Options:
-    region_list: each line can be a region string or a region file
-
-Version: $version
-);
-    print $text;
-}
-## usage }}}
-
 # Example fit config {{{
 #----------------------------------------------------------------
 #    abund grsa
@@ -146,82 +129,91 @@ my $grouptypeval = "20";
 my $binspec = "1:128:2,129:256:4,257:512:8,513:1024:16";
 ## variables }}}
 
+
+## usage {{{
+sub usage {
+    my $prog = basename($0);
+    my $text = qq(SFR - Spectral Fitting Robot
+Usage:
+    $prog <region_list> <evt> <bkgd> <asol> <msk> <bpix> <output> <fit_config>
+
+Options:
+    region_list: each line can be a region string or a region file
+    evt: input event file (support multiple event files with comma ','
+         separated)
+    bkgd: input background file (support blanksky FITS file or
+          local background spectrum file; same number as 'evt')
+    asol: apsect solution files (e.g., '\@asol.lis'; same number as 'evt')
+    msk: detector mask file (same number as 'evt')
+    bpix: bad pixel file (same number as 'evt')
+    output: output text file to store fitting results
+    fit_config: config file that specifies model and output parameters
+
+Version: $version
+);
+    print $text;
+}
+## usage }}}
+
 # 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) {
+my $argc = @ARGV;
+if ($argc != 8) {
     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];
+my $region_list_fn  = $ARGV[0];
+my $evt_fn          = $ARGV[1];
+my $bkg_fn          = $ARGV[2];
+my $asol_fn         = $ARGV[3];
+my $msk_fn          = $ARGV[4];
+my $bpix_fn         = $ARGV[5];
+my $output_fn       = $ARGV[6];
+my $config_fn       = $ARGV[7];
 
 print "====================================================\n";
 print "region_list: $region_list_fn\n";
-print "evt2:        $evt2_fn\n";
-print "bkg:         $bkg_fn \n";
+print "evt:         $evt_fn\n";
+print "bkg:         $bkg_fn\n";
+print "asol:        $asol_fn\n";
+print "msk:         $msk_fn\n";
+print "bpix:        $bpix_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;
+my @evt_list  = split(/[,\ ]+/, $evt_fn);
+my $evt_num   = @evt_list;
+my @bkg_list  = split(/[,\ ]+/, $bkg_fn);
+my $bkg_num   = @bkg_list;
+my @asol_list = split(/[,\ ]+/, $asol_fn);
+my $asol_num  = @asol_list;
+my @msk_list  = split(/[,\ ]+/, $msk_fn);
+my $msk_num   = @msk_list;
+my @bpix_list = split(/[,\ ]+/, $bpix_fn);
+my $bpix_num  = @bpix_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";
+print "evt_num:  $evt_num\n";
+print "bkg_num:  $bkg_num\n";
+print "asol_num: $asol_num\n";
+print "msk_num:  $msk_num\n";
+print "bpix_num: $bpix_num\n";
+if (($evt_num  != $bkg_num) || ($bkg_num != $asol_num) ||
+    ($asol_num != $msk_num) || ($msk_num != $bpix_num)) {
+    print "ERROR: numbers of evt, bkg, asol, msk, bpix 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];
+    printf "--- Data group %d: ---\n", ($i+1);
+    printf "evt:  %s\n", $evt_list[$i];
+    printf "bkg:  %s\n", $bkg_list[$i];
+    printf "asol: %s\n", $asol_list[$i];
+    printf "msk:  %s\n", $msk_list[$i];
+    printf "bpix: %s\n", $bpix_list[$i];
 }
 print "====================================================\n";
-# parse evt, bkg ... }}}
 
 # Backup the output file if it exists
 if ( -e $output_fn ) {
@@ -229,6 +221,25 @@ if ( -e $output_fn ) {
     system("mv -fv $output_fn ${output_fn}_bak");
 }
 
+## 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 .");
+}
+
+$ENV{PFILES} = "./:$ENV{PFILES}";
+## PFILES }}}
+
 # Determine the type of each background file
 my @bkg_type = ();
 foreach my $bkgd (@bkg_list) {
@@ -254,14 +265,20 @@ if (-e $fit_result_fn) {
     unlink($fit_result_fn);
 }
 
+my @fit_outroot_list = ();
 my @fit_spec_list = ();
 my @fit_grp_spec_list = ();
 my @fit_bkg_spec_list = ();
+my @fit_arf_list = ();
+my @fit_rmf_list = ();
 foreach my $i (0 .. $#evt_list) {
     my $ii = $i + 1;
+    push (@fit_outroot_list,  "_fit_${ii}");
     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");
+    push (@fit_arf_list,      "_fit_${ii}.arf");
+    push (@fit_rmf_list,      "_fit_${ii}.rmf");
 }
 ## varaibles }}}
 
@@ -330,7 +347,6 @@ sub bkg_rescale {
     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`;
@@ -366,13 +382,49 @@ sub bkg_rescale {
 }
 # rescale }}}
 
+# Extract spectrum and make ARF & RMF {{{
+# Usage: make_spectrum($input_fits, $region_file, $outroot, $asp, $msk, $bpix)
+sub make_spectrum {
+    my $argc = @_;
+    if ($argc != 6) {
+        print "ERROR: sub make_spectrum() requires 6 arguments!\n";
+        exit 31;
+    }
+    my $infile   = $_[0];
+    my $regfile  = $_[1];
+    my $outroot  = $_[2];
+    my $asol_lis = $_[3];
+    my $mskfile  = $_[4];
+    my $bpixfile = $_[5];
+    my $cmd_str = "punlearn specextract && "
+        . "specextract infile=\"$infile\[sky=region\($regfile\)\]\" "
+        . "outroot=\"$outroot\" bkgfile=\"\" "
+        . "asp=\"$asol_lis\" mskfile=\"$mskfile\" "
+        . "badpixfile=\"$bpixfile\" "
+        . "correctpsf=no weight=yes weight_rmf=yes "
+        . "energy=0.3:11.0:0.01 channel=1:1024:1 "
+        . "bkgresp=no combine=no binarfwmap=2 "
+        . "grouptype=NONE binspec=NONE "
+        . "verbose=1 clobber=yes";
+    #print $cmd_str;
+    system($cmd_str);
+
+    if (! -e "$outroot.rmf") {
+        symlink("$outroot.wrmf", "$outroot.rmf");
+    }
+    if (! -e "$outroot.arf") {
+        symlink("$outroot.warf", "$outroot.arf");
+    }
+}
+# make spectrum }}}
+
 # 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;
+        print "ERROR: sub extract_spectrum() requires 3 arguments!\n";
+        exit 32;
     }
     my $infile  = $_[0];
     my $outfile = $_[1];
@@ -415,9 +467,9 @@ sub parse_result ($\@\%) {
         print "ERROR: sub parse_result() requires 3 argument!\n";
         exit 51;
     }
-    my $line=$_[0];
+    my $line = $_[0];
     my @para_name_list = @{$_[1]};
-    my %para_num= %{$_[2]};
+    my %para_num = %{$_[2]};
 
     my @parsed_results = ();
 
@@ -457,8 +509,8 @@ open(my $FIT_SCRIPT, ">", $fit_script_fn)
 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 "response 1:${ii} $fit_rmf_list[$i]\n";
+    print $FIT_SCRIPT "arf 1:${ii} $fit_arf_list[$i]\n";
     print $FIT_SCRIPT "backgrnd ${ii} $fit_bkg_spec_list[$i]\n";
 }
 
@@ -643,10 +695,12 @@ foreach my $line (@region_lines) {
     # 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);
+        # Prepare source spectrum, ARF, and RMF
+        make_spectrum($evt_list[$i], $fit_reg_fn,
+                      $fit_outroot_list[$i], $asol_list[$i],
+                      $msk_list[$i], $bpix_list[$i]);
         group_spectrum($fit_spec_list[$i], $fit_grp_spec_list[$i],
-            $grouptype, $grouptypeval, $binspec);
+                       $grouptype, $grouptypeval, $binspec);
         # Prepare background spectrum
         if ($bkg_type[$i] == NO_BKG) {
             # Dose not use a background
@@ -663,7 +717,7 @@ foreach my $line (@region_lines) {
             }
         } elsif ($bkg_type[$i] == BKG_SPEC) {
             # Directly use given background spectrum
-            system("cp -v $bkg_list[$i] $fit_bkg_spec_list[$i]");
+            copy($bkg_list[$i], $fit_bkg_spec_list[$i]);
             if ($rescale) {
                 bkg_rescale($fit_spec_list[$i], $fit_bkg_spec_list[$i]);
             }