/Xspecscript_tie_bin.xcm

  1log xspec_temp.log
  2set xs_return_result 1
  3# Fit X-ray spectral annuli in X-spec
  4# additional functionality to tie parameters together such as metallicity in annuli pairs to give the same fit
  5
  6
  7#Filename inputs
  8
  9#####Check for input file
 10##put "Enter filename with comma separated data obsid, etc."
 11##gets stdin filename
 12
 13
 14# Inputs:
 15# List of Spectra - currently constructs based on prefix and suffix 
 16# Redshift
 17# nH in 10^22 cm^2
 18# Tx guess
 19# Fe guess
 20# * combine bin and unbinned and set stat and binning parameters
 21
 22
 23#set check 0  # if there was a typo doesn't crash
 24#while {$check<1} {
 25#puts "Enter Obs ID: "
 26#gets stdin obsid
 27#puts "Enter number of spectra to fit: "
 28#gets stdin numspec
 29#puts "Enter file prefix (e.g. ann): "
 30#gets stdin prefix
 31#puts "Enter file suffix (e.g. .ai.fits): "
 32#gets stdin suffix
 33##############Enter parameters
 34#puts "Enter redshift: "
 35#gets stdin red 
 36#puts "Enter nH column in units of 10^22 cm^2: "
 37#gets stdin nH
 38#puts "Enter Tx guess: "
 39#gets stdin Tx
 40#puts "Enter Abundance guess: "
 41#gets stdin Fe
 42#puts "The obsid is $obsid"
 43#puts "The number of spectra to be fit is $numspec"
 44#puts "The first file to be fit is $prefix\_$obsid\_0$suffix"
 45#puts "If this first file looks correct enter 1, otherwise enter 0"
 46#gets stdin check
 47#} 
 48
 49#################
 50
 51#Fit Xspec data
 52statistic cstat
 53weight standard
 54method simplex
 55setplot energy
 56setplot xlog on
 57setplot ylog on
 58# Keep going until fit converges.
 59query yes
 60# Open the file to put the results in.
 61set fileid [open fit_result.dat w]
 62cosmo 70 0 0.7
 63#Format is nH, kT, nH(fixed at 1.0), Abundance, redshift, switch, and normalize
 64abund grsa
 65
 66set prefix ann
 67set obsid 11724
 68set suffix .pi
 69set nH 0.010
 70set Tx 7.1
 71set Fe 0.3
 72set red 0.213
 73set Tsoft 0.4
 74set Fesoft 0.3
 75set numspec 9
 76
 77set groups "1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40"
 78#load all the annuli first so regions can be tied across
 79set y 1
 80setplot device $prefix\_$obsid\_bin-all.ps/ps
 81for {set x 0} {$x < $numspec} { incr x } {
 82	data $y:$y $prefix\_$obsid\_$x$suffix	
 83	set curr_group [string range $groups 0 [expr {2 * $x}]]
 84	puts "$curr_group"
 85	setplot group $curr_group
 86	plot data
 87	incr y
 88}
 89
 90# focus on 0.7-7.0 keV range only need to constrain spectra, not groups	
 91ignore bad	
 92ignore 1-$numspec:**-0.7 7.0-** 
 93
 94# input initial values for all the spectra phabs(mekal_Tx+mekal_softXraybg)
 95# this setup automatically ties all parameters across spectra
 96# total of 13 parameters
 97model phabs(mekal+mekal) & $nH & $Tx & 1.0 & $Fe,0 & $red & 0 & 0.01 & $Tsoft & 1.0 & $Fesoft,0 & $red & 0 & 0.01
 98# freeze: nH phabs, Tsoft. Rest are already frozen
 99freeze 1 
100freeze 8
101freeze 10
102#untie: all Tx, all norm, central Fe, pair Fe after that
103set paramlen 13
104set Fe1 17
105set Fe2 30
106set norm1 20
107set norm2 26
108set Tx1 15
109for {set x 0} {$x < [expr {$numspec - 1}]} { incr x } {
110	untie $Tx1 $norm1 $norm2
111	set norm1 [expr {$norm1 + $paramlen}]
112	set norm2 [expr {$norm2 + $paramlen}]
113	set Tx1 [expr {$Tx1 + $paramlen}]
114	
115	if {[expr {$x % 2}] == 0} {
116		newpar $Fe2=$Fe1
117		set Fe1 [expr {$Fe1 + $paramlen}]
118		set Fe2 [expr {$Fe2 + $paramlen}]
119	}
120	
121}
122puts "Waiting here before fit. Press return to continue."
123gets stdin empty
124
125# fit the spectra
126fit 10000 0.001
127fit
128# compute errors on all free parameters	
129# output data to files
130
131for {set x 0} {$x < $numspec} { incr x } {
132	error [expr {1 + $x * $paramlen}]-[expr {13 + $x * $paramlen}]
133	puts "#### Hey Spectra # $x"
134	tclout param [expr {2 + $x * $paramlen}]
135	puts "#### Hey tx $xspec_tclout"
136	tclout param [expr {4 + $x * $paramlen}]
137	puts "#### Hey fe $xspec_tclout"
138	tclout param [expr {5 + $x * $paramlen}]
139	puts "#### Hey z $xspec_tclout"
140	tclout stat
141	puts "#### Hey stat $xspec_tclout"
142	tclout dof
143	puts "Hey dof $xspec_tclout"
144	setplot energy
145	setplot xlog on
146	setplot ylog on
147	setplot 
148	setplot device $prefix\_$obsid\_$x\_bin.ps/ps
149	plot data resid
150}
151close $fileid
152log none
153file copy -force xspec_temp.log xspec_tx_fit_$obsid\_bin.log
154file delete xspec_temp.log
155exit