From b82305558b8c574092aaf6bfda515190918f85bc Mon Sep 17 00:00:00 2001
From: Aaron LI <aly@aaronly.me>
Date: Sat, 18 Nov 2017 11:44:42 +0800
Subject: Add xspec_instlines_mc.tcl to evaluate instrumental lines systematic
effects
See the script header for more explanations.
---
astro/xmm/xspec_instlines_mc.tcl | 281 +++++++++++++++++++++++++++++++++++++++
1 file changed, 281 insertions(+)
create mode 100644 astro/xmm/xspec_instlines_mc.tcl
diff --git a/astro/xmm/xspec_instlines_mc.tcl b/astro/xmm/xspec_instlines_mc.tcl
new file mode 100644
index 0000000..4a3d579
--- /dev/null
+++ b/astro/xmm/xspec_instlines_mc.tcl
@@ -0,0 +1,281 @@
+#
+# Copyright (c) 2017 Weitian LI <weitian@aaronly.me>
+# MIT License
+#
+# XSPEC (v12) TCL script
+#
+# Description
+# -----------
+# This script employs the Monte Carlo method to evaluate the impacts of
+# the instrumental lines (modeled as Gaussians) on the fitting results.
+# For each iteration, the norms of the Gaussians are randomly altered
+# according to the its fitted value and sigma, and fixed, and then fit
+# the model to obtain the new result. All Monte Carlo results are saved
+# to a text file for later analyses, in order to evaluate the systematic
+# impacts of the instrumental lines.
+#
+# Usage
+# -----
+# 1. already fitted spectral model (proper parameters free/thaw/link etc.)
+# 2. "set mc_times <mc-times>"; number of Monte Carlo times (default: 1000)
+# 3. "set outfile <outfile>"; output file to save the results
+# 4. @<this-script>
+#
+# References
+# ----------
+# * XSPEC - tclout, tcloutr
+# https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XStclout.html
+# https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XStcloutr.html
+# * XSPEC - The User Interface
+# https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XSappendixTcl.html
+# * TCL 8.5 - string
+# https://www.tcl.tk/man/tcl8.5/TclCmd/string.htm
+#
+# 2017-11-17
+#
+
+#
+# XSPEC control
+#
+
+# Do not echo commands
+set xs_echo_script 0
+# Renormalize only at the beginning of a fit
+renorm prefit
+# Do not ask whether to continue, just continue
+set query [ tcloutr query ]
+query yes
+# Do not chatter too much
+set chatlevel [ lindex [ tcloutr chatter ] 0 ]
+chatter 5
+
+#
+# Settings
+#
+
+# Name of this script
+set script_name "xspec_instlines_mc.tcl"
+set date_now [ clock format [ clock seconds ] -format %Y%m%d ]
+# Default output filename to store the Monte Carlo results
+set outfile_default "xspec_instlines_mc.${date_now}.csv"
+# Default file to log the XSPEC messages (as the screen output is suppressed)
+set logfile_default "xspec_instlines_mc.${date_now}.log"
+# Default number of Monte Carlo times
+set mc_times_default 1000
+
+#
+# Procedures
+#
+
+# Save current XSPEC fitting results
+proc save_xspec {} {
+ set now [ clock format [ clock seconds ] -format %Y%m%d%H%M ]
+ set xspec_outfile "xspec_saveall.${now}.xcm"
+ save all $xspec_outfile
+}
+
+# Print the header line of the given parameters
+proc print_header {parameters {fd stdout} {sep ,}} {
+ set plabels {}
+ foreach p $parameters {
+ set dg_ [ dict get $p datagrp ]
+ set cname_ [ dict get $p comp_name ]
+ set cnum_ [ dict get $p comp_num ]
+ set pname_ [ dict get $p par_name ]
+ set pnum_ [ dict get $p par_num ]
+ set label "${dg_}/${cname_}.${cnum_}/${pname_}(${pnum_})"
+ lappend plabels $label
+ }
+ puts $fd [ join $plabels $sep ]
+}
+
+# Print the data line of the given parameters
+proc print_data {parameters {fd stdout} {sep ,}} {
+ set pvalues {}
+ foreach p $parameters {
+ set pnum_ [ dict get $p par_num ]
+ set pval_ [ lindex [ tcloutr param $pnum_ ] 0 ]
+ lappend pvalues $pval_
+ }
+ puts $fd [ join $pvalues $sep ]
+}
+
+# Generate normal distribution random number
+# Box-Mueller method
+proc prng_normal {{mean 0.0} {stdev 1.0}} {
+ set pi 3.141592653589793
+ set rad_ [ expr {sqrt(-2.8 * log(rand()))} ]
+ set phi_ [ expr {2.0 * $pi * rand()} ]
+ set r_ [ expr {$rad_ * cos($phi_)} ]
+ return [ expr {$mean + $stdev * $r_} ]
+}
+
+# Randomize the norms of the Gaussian models according to their fitted
+# values and sigmas
+proc randomize_norms {parameters} {
+ foreach p $parameters {
+ set pnum_ [ dict get $p norm_par ]
+ set value_ [ dict get $p norm_value ]
+ set sigma_ [ dict get $p norm_sigma ]
+ if {$sigma_ > 0.0} {
+ set newval_ [ prng_normal $value_ $sigma_ ]
+ if {$newval_ < 0.0} {
+ set newval_ 0.0
+ }
+ newpar $pnum_ $newval_
+ }
+ }
+}
+
+# Freeze the norms of the Gaussians models
+proc freeze_norms {parameters} {
+ foreach p $parameters {
+ set pnum_ [ dict get $p norm_par ]
+ set pdelta_ [ lindex [ tcloutr param $pnum_ ] 1 ]
+ set plink_ [ lindex [ tcloutr plink $pnum_ ] 0 ]
+ if {$pdelta_ > 0.0 && $plink_ == "F"} {
+ freeze $pnum_
+ }
+ }
+}
+
+#
+# Main
+#
+
+save_xspec
+
+# Output log file
+if {[ info exists logfile ]} {
+ puts "Output results file: ${logfile}"
+} else {
+ puts "WARNING: logfile not set, default to ${logfile_default}"
+ set logfile $logfile_default
+}
+if {[ file exists $logfile ]} {
+ file rename -force $logfile ${logfile}.old
+}
+# Enable log to file
+log $logfile
+
+# Output results file
+if {[ info exists outfile ]} {
+ puts "Output results file: ${outfile}"
+} else {
+ puts "WARNING: outfile not set, default to ${outfile_default}"
+ set outfile $outfile_default
+}
+if {[ file exists $outfile ]} {
+ file rename -force $outfile ${outfile}.old
+}
+set outfd [ open $outfile w ]
+
+if {[ info exists mc_times ]} {
+ puts "Number of Monte Carlo iterations: ${mc_times}"
+} else {
+ puts "WARNING: mc_times not set, default to ${mc_times_default} times"
+ set mc_times $mc_times_default
+}
+
+set ndatagrp [ tcloutr datagrp ]
+set nmodcomp [ tcloutr modcomp ]
+puts "Number of data groups: ${ndatagrp}"
+puts "Number of model components: ${nmodcomp}"
+
+# List of Gaussian lines, of which each element is an dictionary of elements:
+# - datagrp : data group number
+# - comp : model component number
+# - norm_par : parameter number of the Gaussian norm
+# - norm_value : Gaussian norm value
+# - norm_sigma : Gaussian norm sigma (i.e., standard deviation)
+set gaussians {}
+puts "-----------------------------------------------------------------------"
+for {set i 1} {$i <= $ndatagrp} {incr i} {
+ for {set j 1} {$j <= $nmodcomp} {incr j} {
+ set comp_ [ tcloutr compinfo $j $i ]
+ set comp_name_ [ lindex $comp_ 0 ]
+ if {[ string equal $comp_name_ gaussian ]} {
+ set norm_par_ [ expr {[ lindex $comp_ 1 ] + 2}]
+ set norm_value_ [ lindex [ tcloutr param $norm_par_ ] 0 ]
+ set norm_sigma_ [ tcloutr sigma $norm_par_ ]
+ set gaus_ [ dict create \
+ datagrp $i \
+ comp $j \
+ norm_par $norm_par_ \
+ norm_value $norm_value_ \
+ norm_sigma $norm_sigma_ ]
+ lappend gaussians $gaus_
+ puts "${i}/gaussian: ${gaus_}"
+ }
+ }
+}
+puts "Number of Gaussian lines: [ llength $gaussians ]"
+
+# List of free parameters, of which each element is an dictionary of elements:
+# - datagrp : data group number
+# - comp_num : model component number
+# - comp_name : model component name
+# - par_num : parameter number
+# - par_name : parameter name
+set freeparameters {}
+puts "-----------------------------------------------------------------------"
+for {set i 1} {$i <= $ndatagrp} {incr i} {
+ for {set j 1} {$j <= $nmodcomp} {incr j} {
+ set comp_ [ tcloutr compinfo $j $i ]
+ set comp_name_ [ lindex $comp_ 0 ]
+ set comp_pfirst_ [ lindex $comp_ 1 ] ;# 1st parameter number
+ set comp_np_ [ lindex $comp_ 2 ] ;# number of parameters
+ set comp_plast_ [ expr {$comp_pfirst_ + $comp_np_ - 1} ]
+ for {set k $comp_pfirst_} {$k <= $comp_plast_} {incr k} {
+ scan [ tcloutr param $k ] "%f %f" par_value_ par_delta_
+ set par_link_ [ lindex [ tcloutr plink $k ] 0 ]
+ set par_sigma_ [ tcloutr sigma $k ]
+ if {$par_delta_ > 0.0 && $par_link_ == "F"} {
+ set par_name_ [ lindex [ tcloutr pinfo $k ] 0 ]
+ set param_ [ dict create \
+ datagrp $i \
+ comp_num $j \
+ comp_name $comp_name_ \
+ par_num $k \
+ par_name $par_name_ ]
+ lappend freeparameters $param_
+ puts "${i}/${comp_name_}.${j}/${par_name_}(${k}):\
+ ${par_value_} +/- ${par_sigma_}"
+ }
+ }
+ }
+}
+puts "Number of free parameters: [ llength $freeparameters ]"
+print_header $freeparameters $outfd
+
+puts "-----------------------------------------------------------------------"
+set tstart [ clock seconds ]
+for {set i 0} {$i < $mc_times} {incr i} {
+ puts -nonewline "... [ expr {$i + 1} ] / ${mc_times} ..."
+ if {$i == 0} {
+ puts ""
+ } else {
+ set tnow [ clock seconds ]
+ set elapsed [ expr {($tnow - $tstart) / 60.0} ]
+ set eta [ expr {$elapsed * ($mc_times - $i) / $i} ]
+ puts [ format " Elapsed %.1f min / ETA %.1f min ..." $elapsed $eta ]
+ }
+ randomize_norms $gaussians
+ freeze_norms $gaussians
+ fit
+ print_data $freeparameters $outfd
+}
+puts "-----------------------------------------------------------------------"
+set tnow [ clock seconds ]
+set elapsed [ expr {($tnow - $tstart) / 60.0} ]
+puts [ format "Total elapsed time: %.1f min" $elapsed ]
+
+close $outfd
+# Recover query and chatter level
+query $query
+chatter $chatlevel
+# Disable log to file
+log none
+
+puts "Check the log file for more details: ${logfile}"
+puts "Results written into: ${outfile}"
--
cgit v1.2.2