PhotonStatsAnalysis/StandardDeviation.py at main · JCodeZero/PhotonStatsAnalysis · 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
import sys
from astropy.io import fits
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

filename = sys.argv[1]
eventfile=fits.open(filename)
events=eventfile[1].data
shortdat=events[(events['TIME']>660920900) & (events['TIME']<660921200)]
cleanevts = shortdat[(shortdat['EVENT_FLAGS']==0) & (shortdat['ENERGY']>15)]
filterdat = shortdat[(shortdat['EVENT_FLAGS']==0) & (shortdat['ENERGY']>15)]

time = []

df = pd.DataFrame(filterdat, columns = ['TIME','DET_ID','EVENT_FLAGS','PHA','MASK_WEIGHT','DETX','DETY','PI','ENERGY'])
df['TIME'] = df['TIME'].apply(lambda x: round(x))
time.extend(df['TIME'].tolist())

photons = np.empty((0,1), int)
temp = 0

timearray = np.array(list(range(time[0], time[-1])))


for i in timearray:
    for j in range(len(time)):
        if time[j] == i:
            temp = temp + 1
    photons = np.append(photons, np.array([[temp]]))
    temp = 0

length = 5
edgesx = np.empty((0,1))
means = np.empty((0,1))
stds = np.empty((0,1))
topy = np.empty((0,1))
bottomy = np.empty((0,1))
for i in range(int((time[-1] - time[0])/length)):
    edgesx = np.append(edgesx, np.array([[time[0] + length*i]]))
    means = np.append(means, np.array([[np.mean(photons[length*i:length*(i+1)])]]))
    std = np.std(photons[length*i:length*(i+1)])
    stds = np.append(stds, np.array([[std]]))
    topy = np.append(topy, np.array([[np.mean(photons[length*i:length*(i+1)]) + 3 * std]]))
    bottomy = np.append(bottomy, np.array([[np.mean(photons[length*i:length*(i+1)]) - 3 * std]]))

edgesx = np.append(edgesx, np.array([[edgesx[-1] + length]]))
means = np.append(means, np.array([[means[-1]]]))
stds = np.append(stds, np.array([[stds[-1]]]))
topy = np.append(topy, np.array([[topy[-1]]]))
bottomy = np.append(bottomy, np.array([[bottomy[-1]]]))

plt.plot(timearray, photons, 'r-', drawstyle='steps-post')
plt.plot(edgesx, means, 'blue', drawstyle='steps-post')
plt.plot(edgesx, topy, 'green', drawstyle='steps-post')
plt.plot(edgesx, bottomy, 'green', drawstyle='steps-post')
plt.xlim(660920950,660921100)
plt.ylim(0,20000)
plt.show()

plt.plot(timearray, photons, 'r-', drawstyle='steps-post')
plt.plot(edgesx, means, 'blue', drawstyle='steps-post')
plt.plot(edgesx, topy, 'green', drawstyle='steps-post')
plt.plot(edgesx, bottomy, 'green', drawstyle='steps-post')
plt.xlim(660921050,660921200)
plt.ylim(0,20000)
plt.show()