aboutsummaryrefslogtreecommitdiffstats
path: root/astro/spectrum
diff options
context:
space:
mode:
authorAaron LI <aly@aaronly.me>2019-10-21 22:29:34 +0800
committerAaron LI <aly@aaronly.me>2019-10-21 22:29:34 +0800
commit70a7ec1719bf5ee9cedac46412fb4f7dbb9de926 (patch)
tree4d4600c3dadbcefe344d65952dc440f32fe70138 /astro/spectrum
parent7a626f6e242adf473184d4c0bb06929e5ce2fb43 (diff)
downloadatoolbox-70a7ec1719bf5ee9cedac46412fb4f7dbb9de926.tar.bz2
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
Diffstat (limited to 'astro/spectrum')
-rwxr-xr-xastro/spectrum/sfr_chandra.pl232
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]);
}