astro/xmm: Add several previously developed scripts

1 files changed, 447 insertions, 0 deletions

diff --git a/astro/xmm/xspec_bkgcorr_xmm.tcl b/astro/xmm/xspec_bkgcorr_xmm.tcl
new file mode 100644
index 0000000..240d4a2
--- /dev/null
+++ b/astro/xmm/xspec_bkgcorr_xmm.tcl
@@ -0,0 +1,447 @@
+#####################################################################
+## XSPEC TCL script
+##
+## To apply background spectrum correction to MOS1 and MOS2 background.
+## It fakes the XB (Galatic + CXB) and NXB (SP + instrument lines)
+## spectra according to the current fitting results, and then invoke
+## "mathpha" to make the final corrected background spectrum.
+##
+## NOTE:
+## The exposure time for the faked spectra is multiplied by a scale
+## factor in order to obtain satisfactory statistics for the higher
+## energy regime.  The same scale factor is also applied to the original
+## FWC background spectrum, considering that many channels have "counts"
+## smaller than one, especially in the higher energy regime.  Otherwise,
+## many of these values will be rounded to *zero* by the "mathpha"
+## operating in the "COUNT" mode.
+## See also the following description to variable "SCALE_FACTOR",
+## procedures "spectrum_rate2counts" and "scale_spectrum".
+##
+##
+## NOTES:
+## Work with XSPEC v12.x
+## Based on "xspec_bkgcorr_v2.tcl" for Chandra background spectrum correction
+##
+## Weitian LI <liweitianux@gmail.com>
+## Created: 2015-11-07
+## Updated: 2015-11-12
+##
+## ChangeLogs:
+## 2015-11-10:
+##   * Remove "diag_rsp"
+##   * Add proc update_header:
+##     Update specified keywords instead of just copy the original header
+##   * Use lists to record the information for all data groups,
+##     therefore this tool can handle any number of data groups.
+## 2015-11-11:
+##   * Add back "diag_rsp"
+## 2015-11-12:
+##   * Update to use "ftkeypar" and "fthedit"
+##   * Add procedure "scale_spectrum" to scale the FWC background
+##   * Add procedure "make_filename"
+##   * Add procedure "spectrum_rate2counts" to convert RATE spectrum to
+##     COUNTS, therefore this script can properly handle PN spectrum
+##   * Update script description
+#####################################################################
+
+
+## Global variables {{{
+set datagrp_total [ tcloutr datagrp ]
+
+# Lists to record the information of each data group
+set rootname       {}
+set rmf            {}
+set arf            {}
+set diag_rsp       {}
+set src_spec       {}
+set bkg_spec       {}
+set src_exposure   {}
+set bkg_exposure   {}
+set src_backscale  {}
+set bkg_backscale  {}
+set modelpars      {}
+set back_modelname {}
+set back_modelcomp {}
+set back_modelpars {}
+# This "modelcomp" is not a list
+set modelcomp      {}
+
+# FITS header keywords needed to be copyed to the corrected spectrum
+set KEYWORDS       {TELESCOP INSTRUME FILTER}
+
+# Scale factor for faking spectra
+# By lengthen the exposure time for the fake process, the total counts
+# of the faked spectra is much larger, therefore the statistical issuses
+# due to small counts (especially in the higher energy regime) is mitigated.
+set SCALE_FACTOR   100
+## Global variables }}}
+
+## Procedures
+# Make new filename by appending the suffix to the stem of filename.
+# e.g., "stem.ext" + "sfix" => "stem_sfix.ext"
+proc make_filename {filename suffix} {
+    regexp -- {(.+)\.([^\.]+)} $filename match stem_ ext_
+    set newname "${stem_}_${suffix}.${ext_}"
+    return $newname
+}
+
+# Determine root name and get other filenames (e.g., rmf, arf, back)
+proc get_filenames {datagrp} {
+    global rootname
+    global src_spec
+    global bkg_spec
+    global rmf
+    global arf
+    global diag_rsp
+    set src_spec_ [ exec basename [ tcloutr filename $datagrp ] ]
+    regexp -- {(.+)_(bin|grp|group)([0-9]*)\.(pi|fits)} $src_spec_ match rootname_ grptype_ grpval_ fext_
+    set bkg_spec_ [ tcloutr backgrnd $datagrp ]
+    set response_ [ tcloutr response $datagrp ]
+    set rmf_      [ lindex $response_ 0 ]
+    set diag_rsp_ [ lindex $response_ 1 ]
+    set arf_      [ lindex [ tcloutr arf      $datagrp ] 0 ]
+    lappend rootname $rootname_
+    lappend src_spec $src_spec_
+    lappend bkg_spec $bkg_spec_
+    lappend rmf      $rmf_
+    lappend arf      $arf_
+    lappend diag_rsp $diag_rsp_
+}
+
+# Get exposure and backscale values
+proc get_exposure_backscale {datagrp} {
+    global src_exposure
+    global bkg_exposure
+    global src_backscale
+    global bkg_backscale
+    set src_exposure_ [ tcloutr expos $datagrp s ]
+    set bkg_exposure_ [ tcloutr expos $datagrp b ]
+    scan [ tcloutr backscal $datagrp s ] "%f" src_backscale_
+    scan [ tcloutr backscal $datagrp b ] "%f" bkg_backscale_
+    lappend src_exposure  $src_exposure_
+    lappend bkg_exposure  $bkg_exposure_
+    lappend src_backscale $src_backscale_
+    lappend bkg_backscale $bkg_backscale_
+}
+
+# Get the model components
+proc get_modelcomp {} {
+    global datagrp_total
+    global modelcomp
+    # NOTE: "model" contains all the model components corresponding to
+    # echo data group, and contains also the background model components.
+    # Therefore, the "model" contents repeats the same model components
+    # for each data group.
+    # e.g.,
+    # if data group 1 has model components: "wabs*apec", and there are
+    # altogether 2 data groups, then the "model" content is:
+    #     "wabs(apec) wabs(apec)"
+    # which is repeated!
+    set model [ tcloutr model ]
+    set indices [ lsearch -all $model {*:} ]
+    if { [ llength $indices ] == 0 } {
+        set modelcomp $model
+    } else {
+        set modelcomp [ lrange $model 0 [ expr {[ lindex $indices 0 ] - 1} ] ]
+    }
+    # Deal with the repeated model components if exist multiple data groups
+    set comp_length [ expr { [ llength $modelcomp ] / $datagrp_total } ]
+    set modelcomp [ lrange $modelcomp 0 [ expr {$comp_length - 1} ] ]
+}
+
+# Get the model parameters
+proc get_modelpars {datagrp} {
+    global datagrp_total
+    global modelpars
+    set modelpars_ ""
+    # NOTE: "modpar" is the total number of parameters,
+    # not just for one data group
+    set npar_ [ expr { [ tcloutr modpar ] / $datagrp_total } ]
+    set par_begin_ [ expr {$npar_ * ($datagrp - 1) + 1} ]
+    set par_end_   [ expr {$npar_ + $par_begin_ - 1} ]
+    for {set i $par_begin_} {$i <= $par_end_} {incr i} {
+        scan [ tcloutr param $i ] "%f" pval_
+        set modelpars_ "$modelpars_ $pval_ &"
+    }
+    lappend modelpars $modelpars_
+}
+
+# Get the model name, components and parameters of the background model
+# corresponding to the given data group
+proc get_backmodel {datagrp} {
+    global back_modelname
+    global back_modelcomp
+    global back_modelpars
+    set model_ [ tcloutr model ]
+    set indices_ [ lsearch -all $model_ {*:} ]
+    set name_idx_ [ lindex $indices_ [ expr {$datagrp - 1} ] ]
+    set back_modelname_ [ regsub {:\s*} [ lindex $model_ $name_idx_ ] "" ]
+    set ncomp_ [ tcloutr modcomp $back_modelname_ ]
+    set back_modelcomp_ [ lrange $model_ [ expr {$name_idx_ + 1} ] [ expr {$name_idx_ + $ncomp_ } ] ]
+    set npar_ [ tcloutr modpar $back_modelname_ ]
+    set back_modelpars_ ""
+    for {set i 1} {$i <= $npar_} {incr i} {
+        scan [ tcloutr param ${back_modelname_}:$i ] "%f" pval_
+        set back_modelpars_ "$back_modelpars_ $pval_ &"
+    }
+    lappend back_modelname $back_modelname_
+    lappend back_modelcomp $back_modelcomp_
+    lappend back_modelpars $back_modelpars_
+}
+
+# Save current XSPEC fitting results
+proc save_xspec_fit {} {
+    set now [ exec date +%Y%m%d%H%M ]
+    set xspec_outfile "xspec_${now}.xcm"
+    if { [ file exists ${xspec_outfile} ] } {
+        exec mv -fv ${xspec_outfile} ${xspec_outfile}_bak
+    }
+    save all "${xspec_outfile}"
+}
+
+# Load model
+proc load_model {comp pars} {
+    model clear
+    model $comp & $pars /*
+}
+
+# Set the norm of the ICM APEC component to be zero,
+# which is the last parameter of current loaded model.
+proc set_icm_norm {} {
+    set npar_ [ tcloutr modpar ]
+    newpar $npar_ 0.0
+    puts "ICM_apec_norm:$npar_: [ tcloutr param $npar_ ]"
+}
+
+# Fake spectrum
+proc fake_spectrum {outfile rmf arf exptime} {
+    if {[ file exists ${outfile} ]} {
+        exec mv -fv ${outfile} ${outfile}_bak
+    }
+    data none
+    fakeit none & $rmf & $arf & y &  & $outfile & $exptime & /*
+}
+
+# Convert a spectrum of "RATE" to a spectrum of "COUNTS" by multiplying
+# the exposure time.  The "STAT_ERR" column is also scaled if exists.
+# NOTE:
+# The MOS1/2 FWC background spectrum generated by "mos_back" has column
+# "COUNTS", however, the PN FWC background spectrum generated by "pn_back"
+# has column "RATE".
+# Therefore, the PN background spectra should be converted to "COUNTS" format,
+# since the other faked CXB & NXB spectra are also in "COUNTS" format.
+# The exposure time is acquired from the "EXPOSURE" keyword in the header
+# of input spectrum.
+#
+# Return:
+#   * 0  : input spectrum is already in "COUNTS" format, no conversion needed
+#   * 1  : input spectrum is in "RATE" format, and is converted to "COUNTS"
+#          format and write to outfile
+#   * -1 : invalid input spectrum
+proc spectrum_rate2counts {outfile infile} {
+    if {[ file exists ${outfile} ]} {
+        exec mv -fv ${outfile} ${outfile}_bak
+    }
+    # Get number of columns/fields
+    exec ftkeypar "${infile}+1" TFIELDS
+    set fields_ [ exec pget ftkeypar value ]
+    set colnames {}
+    for {set i 1} {$i <= $fields_} {incr i} {
+        exec ftkeypar "${infile}+1" "TTYPE${i}"
+        set val_ [ exec pget ftkeypar value ]
+        set colname_ [ string toupper [ string trim $val_ {' } ] ]
+        lappend colnames $colname_
+    }
+    if { [ lsearch $colnames {COUNTS} ] != -1 } {
+        # Input spectrum is already in "COUNTS" format
+        return 0
+    } elseif { [ lsearch $colnames {RATE} ] != -1 } {
+        # Get exposure time
+        exec ftkeypar "${infile}+1" EXPOSURE
+        set exposure_ [ exec pget ftkeypar value ]
+        # Convert column "RATE"
+        set tmpfile "${outfile}_tmp[pid]"
+        exec cp -fv ${infile} ${tmpfile}
+        # COUNTS = RATE * EXPOSURE
+        set fcmd_ "ftcalc \"${tmpfile}\" \"${outfile}\" column=RATE expression=\"RATE * $exposure_\" history=yes"
+        puts $fcmd_
+        exec sh -c $fcmd_
+        exec rm -fv ${tmpfile}
+        # Update column name from "RATE" to "COUNTS", and column units
+        set field_idx_ [ expr { [ lsearch $colnames {RATE} ] + 1 } ]
+        set fcmd_ "fthedit \"${outfile}+1\" keyword=\"TTYPE${field_idx_}\" operation=add value=\"COUNTS\""
+        puts $fcmd_
+        exec sh -c $fcmd_
+        set fcmd_ "fthedit \"${outfile}+1\" keyword=\"TUNIT${field_idx_}\" operation=add value=\"count\""
+        puts $fcmd_
+        exec sh -c $fcmd_
+        # Convert column "STAT_ERR"
+        if { [ lsearch $colnames {STAT_ERR} ] != -1 } {
+            exec mv -fv ${outfile} ${tmpfile}
+            # STAT_ERR = STAT_ERR * EXPOSURE
+            set fcmd_ "ftcalc \"${tmpfile}\" \"${outfile}\" column=STAT_ERR expression=\"STAT_ERR * $exposure_\" history=yes"
+            puts $fcmd_
+            exec sh -c $fcmd_
+            exec rm -fv ${tmpfile}
+        }
+        return 1
+    } else {
+        # Invalid input spectrum
+        return -1
+    }
+}
+
+# Scale the spectrum "COUNTS" and "STAT_ERR" (if exists) columns by a factor.
+# NOTE:
+# This procedure is mainly used to scale the FWC background, within which
+# many channels have "counts" less than one.  When using "mathpha" to combine
+# this FWC background and other faked spectra in "COUNT" mode, the channels
+# with "counts" less than one (especially the higher energy regime) may be
+# truncated/rounded to integer, therefore the higher energy regime of the
+# combined spectrum may be greatly underestimated, and will cause the
+# following fitting very bad in the higher energy regime.
+# See also the description about the "SCALE_FACTOR" above.
+proc scale_spectrum {outfile infile factor} {
+    if {[ file exists ${outfile} ]} {
+        exec mv -fv ${outfile} ${outfile}_bak
+    }
+    exec cp -fv ${infile} ${outfile}
+    # Get number of columns/fields
+    exec ftkeypar "${outfile}+1" TFIELDS
+    set fields_ [ exec pget ftkeypar value ]
+    for {set i 1} {$i <= $fields_} {incr i} {
+        exec ftkeypar "${outfile}+1" "TTYPE${i}"
+        set val_ [ exec pget ftkeypar value ]
+        set colname_ [ string toupper [ string trim $val_ {' } ] ]
+        # Scale column "COUNTS"
+        if {$colname_ == "COUNTS"} {
+            set tmpfile "${outfile}_tmp[pid]"
+            exec mv -fv ${outfile} ${tmpfile}
+            # COUNTS = COUNTS * factor
+            set fcmd_ "ftcalc \"${tmpfile}\" \"${outfile}\" column=COUNTS expression=\"COUNTS * $factor\" history=yes"
+            puts $fcmd_
+            exec sh -c $fcmd_
+            exec rm -fv ${tmpfile}
+        }
+        # Scale column "STAT_ERR"
+        if {$colname_ == "STAT_ERR"} {
+            set tmpfile "${outfile}_tmp[pid]"
+            exec mv -fv ${outfile} ${tmpfile}
+            # STAT_ERR = STAT_ERR * factor
+            set fcmd_ "ftcalc \"${tmpfile}\" \"${outfile}\" column=STAT_ERR expression=\"STAT_ERR * $factor\" history=yes"
+            puts $fcmd_
+            exec sh -c $fcmd_
+            exec rm -fv ${tmpfile}
+        }
+    }
+}
+
+# Combine faked spectra to original FWC background spectrum with "mathpha"
+proc combine_spectra {outfile fwc_bkg cxb nxb exposure backscale} {
+    set combine_expr_ "$fwc_bkg + $cxb + $nxb"
+    set comment1_ "Corrected background spectrum; FWC based"
+    set combine_cmd_ "mathpha expr='${combine_expr_}' outfil=${outfile} exposure=${exposure} backscal=${backscale} units=C areascal=% properr=yes ncomments=1 comment1='${comment1_}'"
+    if {[ file exists ${outfile} ]} {
+        exec mv -fv ${outfile} ${outfile}_bak
+    }
+    puts $combine_cmd_
+    exec sh -c $combine_cmd_
+}
+
+# Copy header keywords from original background spectrum to the
+# newly combined background spectrum
+proc update_header {newfile origfile} {
+    global KEYWORDS
+    for {set i 0} {$i < [ llength $KEYWORDS ]} {incr i} {
+        set key_ [ lindex $KEYWORDS $i ]
+        exec ftkeypar "${origfile}+1" $key_
+        set val_ [ exec pget ftkeypar value ]
+        #set fcmd_ "fparkey value=\"$val_\" fitsfile=\"${newfile}+1\" keyword=\"$key_\" add=yes"
+        set fcmd_ "fthedit \"${newfile}+1\" keyword=\"$key_\" operation=add value=\"$val_\""
+        puts $fcmd_
+        exec sh -c $fcmd_
+    }
+}
+
+
+# Save current fit results
+save_xspec_fit
+
+# Get the current model components, which is the same for all data groups
+get_modelcomp
+
+for {set dg 1} {$dg <= $datagrp_total} {incr dg} {
+    get_filenames          $dg
+    get_exposure_backscale $dg
+    get_modelpars          $dg
+    get_backmodel          $dg
+}
+
+puts "modelcomp: $modelcomp"
+puts "rootname: $rootname"
+puts "rmf: $rmf"
+puts "arf: $arf"
+puts "diag_rsp: $diag_rsp"
+puts "src_spec: $src_spec"
+puts "bkg_spec: $bkg_spec"
+puts "src_exposure: $src_exposure"
+puts "bkg_exposure: $bkg_exposure"
+puts "src_backscale: $src_backscale"
+puts "bkg_backscale: $bkg_backscale"
+puts "modelpars: $modelpars"
+puts "back_modelname: $back_modelname"
+puts "back_modelcomp: $back_modelcomp"
+puts "back_modelpars: $back_modelpars"
+puts "SCALE_FACTOR: $SCALE_FACTOR"
+
+# DEBUG
+set DEBUG 0
+if {$DEBUG != 1} {
+
+# Clear all loaded model data and models
+data none
+model clear
+
+for {set idg 0} {$idg < $datagrp_total} {incr idg} {
+    set rootname_       [ lindex $rootname       $idg ]
+    set bkg_spec_       [ lindex $bkg_spec       $idg ]
+    set rmf_            [ lindex $rmf            $idg ]
+    set arf_            [ lindex $arf            $idg ]
+    set diag_rsp_       [ lindex $diag_rsp       $idg ]
+    set bkg_exposure_   [ lindex $bkg_exposure   $idg ]
+    set bkg_backscale_  [ lindex $bkg_backscale  $idg ]
+    set modelpars_      [ lindex $modelpars      $idg ]
+    set back_modelcomp_ [ lindex $back_modelcomp $idg ]
+    set back_modelpars_ [ lindex $back_modelpars $idg ]
+    set fake_exposure   [ expr {$bkg_exposure_  * $SCALE_FACTOR} ]
+    set fake_backscale  [ expr {$bkg_backscale_ * $SCALE_FACTOR} ]
+    set fake_cxb        "fake_cxb_${rootname_}.pi"
+    set fake_nxb        "fake_nxb_${rootname_}.pi"
+    set bkgspec_cnts    [ make_filename $bkg_spec_ "cnts" ]
+    set bkgcorr_outfile "bkgcorr_${rootname_}.pi"
+    # Fake CXB and NXB
+    load_model    $modelcomp $modelpars_
+    set_icm_norm
+    fake_spectrum $fake_cxb $rmf_ $arf_ $fake_exposure
+    load_model    $back_modelcomp_ $back_modelpars_
+    fake_spectrum $fake_nxb $diag_rsp_ "" $fake_exposure
+    # Convert FWC background spectrum from "RATE" to "COUNTS" for PN
+    set ret [ spectrum_rate2counts $bkgspec_cnts $bkg_spec_ ]
+    if { $ret == 0 } {
+        set bkgspec_cnts $bkg_spec_
+    } elseif { $ret == 1 } {
+        puts "Converted RATE spectrum to COUNTS: $bkgspec_cnts"
+    } else {
+        return -code error "ERROR: invalid spectrum for 'spectrum_rate2counts'"
+    }
+    set bkgspec_scaled [ make_filename $bkgspec_cnts "sf${SCALE_FACTOR}" ]
+    # Scale original FWC background before combination
+    scale_spectrum  $bkgspec_scaled $bkgspec_cnts $SCALE_FACTOR
+    # Combine faked spectra with original FWC background
+    combine_spectra $bkgcorr_outfile $bkgspec_scaled $fake_cxb $fake_nxb $bkg_exposure_ $fake_backscale
+    update_header   $bkgcorr_outfile $bkg_spec_
+}
+
+}
+
+# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=tcl: #