# t_fits_pl_softf.xcm 4.0
#
# Xspec script
#
# Author: Paul O'Neill
#
# This script fits the spectrum in the energy range 0.6-2.0 keV. The flux is
# calculated in the range 0.5-2.0 keV.
#
# The model DOES NOT CONTAIN an absorption component, and use peg power-law.
#
# The initial fit is shown in the file web_pha_softf_init.gif
#
# The files fit_result_pl_softf.dat, fit_result_flux_softf.dat, web_pha_softf.gif,
# are NOT produced if the Tcl procedure ezyfit_v4.tcl fails.
#
# Calls the script: ezyfit_v4.tcl
#
# Date in final form: 17 December 2004

############################################################################
# Set up and read in data
############################################################################

# Return Tcl results from Xspec commands

set xs_return_result 1

# Set up log file

log plfit.log

# Continue fit without prompting

query no

# Read in the data

data 1:1 SEQUENCE_sis0_src_20
data 2:2 SEQUENCE_sis1_src_20
data 3:3 SEQUENCE_gis2_src_20
data 4:4 SEQUENCE_gis3_src_20

# Ignore bad channels - only use soft band (0.6-2.0 keV)

ignore bad
ignore 1-4:0.1-0.6
ignore 1-4:2.0-20.0

# Set plot parameter

setplot energy

# Source the ezyfit_v4.tcl procedure

source ezyfit_v4.tcl

# Set maximum acceptable chi-squared and number of  fits

set maxaccchi 15.0
set maxnumfit 20


############################################################################
# Define model
############################################################################

model con*pegpwrlw
1.0,-1, 0,0,1,1
2.,0.2,-1,-1,3,3
0.5,-1,0.5,0.5,0.5,0.5
2.0,-1,2.0,2.0,2.0,2.0
10.0,1.0,0,0,1e8,1e8
0.99,0.01,0,0,2,2
=2
=3
=4
=5
0.99,0.01,0,0,2,2
=2
=3
=4
=5
0.99,0.01,0,0,2,2
=2
=3
=4
=5


############################################################################
# Produce a web PHA file showing the initial fit
############################################################################

fit 25

notice 1-4:0.6-10.0
ignore bad

setplot command ma off
setplot command la t "0.6 - 2 keV no absorption - initial fit"
setplot command la file "SEQUENCE web_pha_softf_init"
setplot command time off
setplot command cs 1.3
setplot command h web_pha_softf_init.gif/gif
plot ldata ratio

ignore 1-4:0.1-0.6
ignore 1-4:2.0-20.0

save all SEQUENCE_webfit_softf


############################################################################
# Run ezyfit_v4
############################################################################

# use delta chi-sq = 1 for 68 % confidence intervals 

set fitresult [ezyfit_v4 [list 2 5 6 11 16] [list 1.0 1.0 1.0 1.0 1.0] 0.05 $maxnumfit 25 $maxaccchi 1000]


############################################################################
# Write parameter values to file if ezyfit_v4 was successful
############################################################################

if {[lindex $fitresult 0] == "GOOD"} {

    set chisq [lindex $fitresult 1]
    set dof  [lindex $fitresult 2]
    set redchi [expr $chisq / $dof]

# ezyfit_v4.tcl return string format is (for gamma, for  example): gam lower limit, 
# quality (1=correct or 0=hit limit), best-fitting gamma,
# gam upper limit, quality (1=correct or 0=hit limit)

    set glo [lindex $fitresult 4]
    set gloq [lindex $fitresult 5]
    set gamma [lindex $fitresult 6]
    set ghi [lindex $fitresult 7]
    set ghiq [lindex $fitresult 8]
    set g_err [expr ($ghi - $glo)/2.0]

    set nolo [lindex $fitresult 9]
    set noloq [lindex $fitresult 10]
    set norm [lindex $fitresult 11]
    set nohi [lindex $fitresult 12]
    set nohiq [lindex $fitresult 13]
    set no_err [expr ($nohi - $nolo)/2.0]

    set sis1flo [lindex $fitresult 14]
    set sis1floq [lindex $fitresult 15]
    set sis1f [lindex $fitresult 16]
    set sis1fhi [lindex $fitresult 17]
    set sis1fhiq [lindex $fitresult 18]
    set sis1f_err [expr ($sis1fhi - $sis1flo)/2.0]

    set gis2flo [lindex $fitresult 19]
    set gis2floq [lindex $fitresult 20]
    set gis2f [lindex $fitresult 21]
    set gis2fhi [lindex $fitresult 22]
    set gis2fhiq [lindex $fitresult 23]
    set gis2f_err [expr ($gis2fhi - $gis2flo)/2.0]

    set gis3flo [lindex $fitresult 24]
    set gis3floq [lindex $fitresult 25]
    set gis3f [lindex $fitresult 26]
    set gis3fhi [lindex $fitresult 27]
    set gis3fhiq [lindex $fitresult 28]
    set gis3f_err [expr ($gis3fhi - $gis3flo)/2.0]

    set fileid [open fit_result_pl_softf.dat w+]
    puts $fileid "$chisq $dof $redchi"
    puts $fileid "$glo $gloq $gamma $ghi $ghiq $g_err"
    puts $fileid "$nolo $noloq $norm $nohi $nohiq $no_err"
    puts $fileid "$sis1flo $sis1floq $sis1f $sis1fhi $sis1fhiq $sis1f_err"
    puts $fileid "$gis2flo $gis2floq $gis2f $gis2fhi $gis2fhiq $gis2f_err"
    puts $fileid "$gis3flo $gis3floq $gis3f $gis3fhi $gis3fhiq $gis3f_err"
    close $fileid

}


############################################################################
# Calculate flux and write to file if ezyfit_v4 was successful
############################################################################

# the file fit_result_flux_softf.dat will contain the fluxes of all
# four datasets, on a single line.  The format is:
#
# dataset1_photons/cm^2/s     dataset1_ergs/cm^2/s
#                         ...
# dataset4_photons/cm^2/s     dataset4_ergs/cm^2/s
#
# dummyrsp is used so that the flux can be calculated
# over the 0.5-2 keV range, even if the data don't
# extend down to 0.5 keV.

if {[lindex $fitresult 0] == "GOOD"} {

    dummyrsp 0.5 2.0
    set thefluxes [flux 0.5 2.0]
    response
    set fileid [open fit_result_flux_softf.dat w+]
    puts $fileid "[lindex $thefluxes 0] [lindex $thefluxes 1]"
    puts $fileid "[lindex $thefluxes 2] [lindex $thefluxes 3]"
    puts $fileid "[lindex $thefluxes 4] [lindex $thefluxes 5]"
    puts $fileid "[lindex $thefluxes 6] [lindex $thefluxes 7]"
    close $fileid

}


############################################################################
# Determine background subtracted counting rate if ezyfit_v4 was successful
############################################################################

if {[lindex $fitresult 0] == "GOOD"} {


### get counting rates and error  -  sis0 -> gis3

    tclout rate 1
    set rate1 [lindex $xspec_tclout 0]
    set rate1_e [lindex $xspec_tclout 1]

    tclout rate 2
    set rate2 [lindex $xspec_tclout 0]
    set rate2_e [lindex $xspec_tclout 1]

    tclout rate 3
    set rate3 [lindex $xspec_tclout 0]
    set rate3_e [lindex $xspec_tclout 1]

    tclout rate 4
    set rate4 [lindex $xspec_tclout 0]
    set rate4_e [lindex $xspec_tclout 1]


#### write to file ####

    set fileid [open xspec_rates_softf.dat w+]
    puts $fileid "$rate1 $rate1_e" 
    puts $fileid "$rate2 $rate2_e"
    puts $fileid "$rate3 $rate3_e"
    puts $fileid "$rate4 $rate4_e"
    close $fileid

}


############################################################################
# Produce the web PHA file
############################################################################

if {[lindex $fitresult 0] == "GOOD"} {

    notice 1-4:0.6-10.0
    ignore bad

    show files

    setplot command ma off
    setplot command la t "0.6 - 2 keV no absorption - best fit"
    setplot command la file "SEQUENCE web_pha_softf"
    setplot command time off
    setplot command cs 1.3
    setplot command h web_pha_softf.gif/gif
    plot ldata ratio
}


############################################################################
# All done
############################################################################


exit