diff options
author | Aaron LI <aly@aaronly.me> | 2017-11-18 11:44:42 +0800 |
---|---|---|
committer | Aaron LI <aly@aaronly.me> | 2017-11-18 11:44:42 +0800 |
commit | b82305558b8c574092aaf6bfda515190918f85bc (patch) | |
tree | f7a73246ef3b30d3be39e2fe2927d541fe243e88 /astro/xmm | |
parent | 4037655221275dd8ff1f5f4afa342cf73b68161a (diff) | |
download | atoolbox-b82305558b8c574092aaf6bfda515190918f85bc.tar.bz2 |
Add xspec_instlines_mc.tcl to evaluate instrumental lines systematic effects
See the script header for more explanations.
Diffstat (limited to 'astro/xmm')
-rw-r--r-- | astro/xmm/xspec_instlines_mc.tcl | 281 |
1 files changed, 281 insertions, 0 deletions
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}" |