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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
# Aaron LI
# Created: 2016-01-16
# Updated: 2016-01-16
#
"""
Clean the lightcurve by fitting the RATE data with a Gaussian model,
and discard the time bins with RATE beyond [mean-n*sigma, mean+n*sigma].
"""
__version__ = "0.1.0"
__date__ = "2016-01-16"
import sys
import argparse
from astropy.io import fits
import numpy as np
class LightCurve:
"""
X-ray data light curve class
"""
def __init__(self, lcfile):
f = fits.open(lcfile)
self.lc_data = f[1].data
self.lc_header = f[1].header
self.time = self.lc_data['TIME']
self.rate = self.lc_data['RATE']
self.rate_err = self.lc_data['ERROR']
self.TSTART = self.lc_header['TSTART']
self.TSTOP = self.lc_header['TSTOP']
self.TIMEDEL = self.lc_header['TIMEDEL']
self.TIMEPIXR = self.lc_header['TIMEPIXR']
f.close()
def sigma_clip(self, nsigma=3, maxiter=10):
"""
Iteratively clip the time bins whose value lie beyond the
range [mean-n*sigma, mean+n*sigma].
"""
rate = self.rate
keep_idx = np.ones(rate.shape, dtype=bool) # all True's
keep_num = np.sum(keep_idx)
keep_num0 = np.inf
i = 0
while (keep_num < keep_num0):
if (i >= maxiter):
print("WARNING: maximum iteration limit reached",
file=sys.stderr)
break
keep_num0 = keep_num
i += 1
mean = np.mean(rate[keep_idx])
sigma = np.std(rate[keep_idx])
cut_low = mean - nsigma * sigma
cut_high = mean + nsigma * sigma
keep_idx = np.logical_and((rate >= cut_low), (rate <= cut_high))
keep_num = np.sum(keep_idx)
# save clip results
self.niter = i
self.keep_idx = keep_idx
self.time_clipped = self.time[keep_idx]
self.rate_clipped = self.rate[keep_idx]
def make_gti(self, apply_header=True):
"""
Make new GTIs (good time intervals) according to the clipped
time bins.
"""
frac = 0.01 # TIMEDEL fraction to distingush two time bins
gti_start = []
gti_stop = []
time_start = self.time_clipped
time_stop = time_start + self.TIMEDEL
# first GTI start time
gti_start.append(time_start[0])
for tstart, tstop in zip(time_start[1:], time_stop[:-1]):
if (np.abs(tstart-tstop) <= frac * self.TIMEDEL):
# time bin continues
continue
else:
# a new GTI start
gti_start.append(tstart)
gti_stop.append(tstop)
# last GTI stop time
gti_stop.append(time_stop[-1])
# convert to numpy array
gti_start = np.array(gti_start)
gti_stop = np.array(gti_stop)
if apply_header:
# add TSTART to the time
gti_start += self.TSTART
gti_stop += self.TSTART
# save results
self.gti_start = gti_start
self.gti_stop = gti_stop
def write_gti(self, filename=None, header=True):
"""
Write generated GTIs to file or screen (default)
"""
if isinstance(filename, str):
outfile = open(filename, 'w')
else:
outfile = sys.stdout
#
if header:
outfile.write('# TSTART\tTSTOP\n')
outfile.write('\n'.join([ '%s\t%s' % (tstart, tstop) \
for tstart, tstop in zip(self.gti_start, self.gti_stop) ]))
#
if isinstance(filename, str):
outfile.close()
def main():
parser = argparse.ArgumentParser(
description="Clean light curve by sigma clipping")
parser.add_argument("-V", "--version", action="version",
version="%(prog)s " + "%s (%s)" % (__version__, __date__))
parser.add_argument("infile",
help="input lightcurve file; contains [TIME, RATE] columns")
parser.add_argument("outfile", nargs='?', default=None,
help="output text-format GTI file; for XSELECT filter time")
parser.add_argument("-s", "--nsigma", dest="nsigma", type=float,
default=2.0, help="sigma clipping significant level")
parser.add_argument("-H", "--no-header", dest="noheader",
action="store_true", help="not write header to the output file")
parser.add_argument("-v", "--verbose", dest="verbose",
action="store_true", help="show verbose information")
args = parser.parse_args()
lc = LightCurve(args.infile)
lc.sigma_clip(nsigma=args.nsigma)
lc.make_gti(apply_header=True)
lc.write_gti(filename=args.outfile, header=(not args.noheader))
if args.verbose:
exposure = np.sum(lc.gti_stop - lc.gti_start)
print("# Total GTI: %.2f (s)" % exposure)
if __name__ == "__main__":
main()
# vim: set ts=4 sw=4 tw=0 fenc= ft=python: #
|
