f-test-heasoft-xspec/f-stat.xcm at main · XRByte/f-test-heasoft-xspec · GitHub

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
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
for { set a 1}  {$a < 100001} {incr a} {
	statistic chi
	data 1:1 fake-test-redge-${a}.fak
	backgrnd 1 fake-test-redge-${a}_bkg.fak
	ignore 1:1-3,25-81

	method leven 10 0.01
	abund angr
	xsect vern
	cosmo 70 0 0.73
	xset delta 0.01
	systematic 0.015
	model  TBabs(bbody + redge + nthComp)
              2     -0.001          0          0     100000      1e+06
        1.18204       0.01     0.0001       0.01        100        200
       0.203465       0.01          0          0      1e+20      1e+24
        7.91735      0.001      0.001      0.001        100        100
        1.36916      0.001      0.001      0.001        100        100
      0.0194493       0.01          0          0      1e+20      1e+24
        1.81854       0.01        1.1        1.6          3          4
           1000       -0.1          1          5       1000       1000
= p2
              1         -1          0          0          1          1
              0      -0.01     -0.999     -0.999         10         10
      0.0370075       0.01          0          0      1e+20      1e+24
	query yes
	fit
	tclout stat
	scan $xspec_tclout "%f" chistat_2
	tclout dof
	scan $xspec_tclout "%f" dof_2
	echo $chistat_2
	clear model

	statistic chi
	data 1:1 fake-test-redge-${a}.fak
	backgrnd 1 fake-test-redge-${a}_bkg.fak
	ignore 1:1-3,25-81

	method leven 10 0.01
	abund angr
	xsect vern
	cosmo 70 0 0.73
	xset delta 0.01
	systematic 0.015
	model  TBabs(bbody + nthComp)
		      2     -0.001          0          0     100000      1e+06
		1.23552       0.01     0.0001       0.01        100        200
	       0.199382       0.01          0          0      1e+20      1e+24
		1.59953       0.01        1.1        1.6          3          4
		   1000       -0.1          1          5       1000       1000
	= p2
		      1         -1          0          0          1          1
		      0      -0.01     -0.999     -0.999         10         10
	      0.0204632       0.01          0          0      1e+20      1e+24
	query yes
	fit
	tclout stat
	scan $xspec_tclout "%f" chistat_1
	tclout dof
	scan $xspec_tclout "%f" dof_1
	echo $chistat_1
	model clear

	set num1 [expr $chistat_2 - $chistat_1]
	set den1 [expr $dof_2 - $dof_1]
	set num [expr $num1 * $dof_2]
	set den [expr $den1 * $chistat_2]
	set res [expr $num / $den]
	echo $res >> f_statistic.txt
}