aboutsummaryrefslogtreecommitdiffstats
path: root/astro/xmm/xspec_bkgcorr_xmm.tcl
diff options
context:
space:
mode:
Diffstat (limited to 'astro/xmm/xspec_bkgcorr_xmm.tcl')
-rw-r--r--astro/xmm/xspec_bkgcorr_xmm.tcl447
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: #
Copyright © 2017-2018 Aaron LI. All Rights Reserved.