aboutsummaryrefslogtreecommitdiffstats
path: root/mass_profile/analyze_entropy_profile.py
blob: c1236818a0bef40a733fcf50c60575070f528ebd (plain)
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
#!/usr/bin/env python

import sys
import numpy
import scipy.interpolate

center_entropy_file=open('entropy_center.qdp')
entropy_file=open('summary_entropy.qdp')
confidence_level=.68
rout=float(sys.argv[1])

center_s=0
for i in center_entropy_file:
    r,s=i.split()
    r=float(r)
    s=float(s)
    if r>rout:
        center_s=s
        break

new_data=True


s_list=[]
for i in entropy_file:
    if i[0]=='n':
        new_data=True
        continue
    if new_data==False:
        continue
    r,s=i.split()
    r=float(r)
    s=float(s)
    if r>rout:
        new_data=False
        s_list.append(s)

s_idx=-1

s_list.sort()
for i in range(len(s_list)-1):
    if (center_s-s_list[i])*(center_s-s_list[i+1])<=0:
        m_idx=i
        break


slidx=int(s_idx*(1-confidence_level))
suidx=s_idx-1+int((len(s_list)-s_idx)*confidence_level)


serr1=s_list[slidx]-center_s
serr2=s_list[suidx]-center_s

print("S=\t%e\t %e/+%e keV cm^2 (1 sigma)"%(center_s,serr1,serr2))