1 | #!/usr/bin/python -tt
|
---|
2 | #
|
---|
3 | # Werner Lustermann, Dominik Neise
|
---|
4 | # ETH Zurich, TU Dortmund
|
---|
5 | #
|
---|
6 | import os
|
---|
7 | import sys
|
---|
8 | import numpy as np
|
---|
9 | import ROOT
|
---|
10 | import pylab
|
---|
11 |
|
---|
12 | from scipy.signal import firwin, iirdesign, lfilter
|
---|
13 |
|
---|
14 | ########## BUILDING OF THE SHARED OBJECT FILES ###############################
|
---|
15 | if __name__ == '__main__' and len(sys.argv) > 1 and 'build' in sys.argv[1]:
|
---|
16 | ROOT.gSystem.AddLinkedLibs("-lz")
|
---|
17 | root_make_string = ROOT.gSystem.GetMakeSharedLib()
|
---|
18 | if not "-std=c++0x" in root_make_string:
|
---|
19 | make_string = root_make_string.replace('$Opt', '$Opt -std=c++0x -D HAVE_ZLIB')
|
---|
20 | ROOT.gSystem.SetMakeSharedLib(make_string)
|
---|
21 | ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/izstream.h+O")
|
---|
22 | ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/fits.h+O")
|
---|
23 | ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/zfits.h+O")
|
---|
24 | ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/factfits.h+O")
|
---|
25 | if not "-std=c++0x" in root_make_string:
|
---|
26 | make_string = root_make_string.replace('$Opt', "$Opt -std=c++0x -D HAVE_ZLIB -D'PACKAGE_NAME=\"PACKAGE_NAME\"' "
|
---|
27 | "-D'PACKAGE_VERSION=\"PACKAGE_VERSION\"' -D'REVISION=\"REVISION\"' ")
|
---|
28 | ROOT.gSystem.SetMakeSharedLib(make_string)
|
---|
29 | ROOT.gROOT.ProcessLine(".L extern_Mars_mcore/DrsCalib.h+O")
|
---|
30 |
|
---|
31 | ROOT.gInterpreter.GenerateDictionary("map<string,fits::Entry>","map;string;extern_Mars_mcore/fits.h")
|
---|
32 | ROOT.gInterpreter.GenerateDictionary("pair<string,fits::Entry>","map;string;extern_Mars_mcore/fits.h")
|
---|
33 | ROOT.gInterpreter.GenerateDictionary("map<string,fits::Table::Column>","map;string;extern_Mars_mcore/fits.h")
|
---|
34 | ROOT.gInterpreter.GenerateDictionary("pair<string,fits::Table::Column>","map;string;extern_Mars_mcore/fits.h")
|
---|
35 | ROOT.gInterpreter.GenerateDictionary("vector<DrsCalibrate::Step>","vector;extern_Mars_mcore/DrsCalib.h")
|
---|
36 |
|
---|
37 | ########## USAGE #############################################################
|
---|
38 | if __name__ == '__main__' and len(sys.argv) < 3:
|
---|
39 | print """ Usage:
|
---|
40 | ----------------------------------------------------------------------
|
---|
41 | To just build the shared object libs call:
|
---|
42 | python pyfact.py build
|
---|
43 |
|
---|
44 | To build (if necessary) and open an example file
|
---|
45 | python -i pyfact.py /path/to/data_file.zfits /path/to/calib_file.drs.fits.gz
|
---|
46 |
|
---|
47 | Any of the 3 file 'types': fits.gz zfits fits
|
---|
48 | should be supported.
|
---|
49 |
|
---|
50 | To clean all of automatically built files, do something like:
|
---|
51 | rm *.so *.d AutoDict_* extern_Mars_mcore/*.so extern_Mars_mcore/*.d
|
---|
52 | ----------------------------------------------------------------------
|
---|
53 | """
|
---|
54 | sys.exit(1)
|
---|
55 |
|
---|
56 | ######### START OF PYFACT MODULE ############################################
|
---|
57 | path = os.path.dirname(os.path.realpath(__file__))
|
---|
58 | ROOT.gSystem.Load(path+'/AutoDict_map_string_fits__Entry__cxx.so')
|
---|
59 | ROOT.gSystem.Load(path+'/AutoDict_map_string_fits__Table__Column__cxx.so')
|
---|
60 | ROOT.gSystem.Load(path+'/AutoDict_pair_string_fits__Entry__cxx.so')
|
---|
61 | ROOT.gSystem.Load(path+'/AutoDict_pair_string_fits__Table__Column__cxx.so')
|
---|
62 | ROOT.gSystem.Load(path+'/AutoDict_vector_DrsCalibrate__Step__cxx.so')
|
---|
63 | ROOT.gSystem.Load(path+'/extern_Mars_mcore/fits_h.so')
|
---|
64 | ROOT.gSystem.Load(path+'/extern_Mars_mcore/izstream_h.so')
|
---|
65 | ROOT.gSystem.Load(path+'/extern_Mars_mcore/zfits_h.so')
|
---|
66 | ROOT.gSystem.Load(path+'/extern_Mars_mcore/factfits_h.so')
|
---|
67 | ROOT.gSystem.Load(path+'/extern_Mars_mcore/DrsCalib_h.so')
|
---|
68 | del path
|
---|
69 |
|
---|
70 | class Fits( object ):
|
---|
71 | """ General FITS file access
|
---|
72 |
|
---|
73 | Wrapper for factfits class from Mars/mcore/factfits.h
|
---|
74 | (factfits might not be suited to read any FITS file out there, but certainly
|
---|
75 | all FACT FITS files can be read using this class)
|
---|
76 | """
|
---|
77 | __module__ = 'pyfact'
|
---|
78 | def __init__(self, path):
|
---|
79 | """
|
---|
80 | """
|
---|
81 | if not os.path.exists(path):
|
---|
82 | raise IOError(path+' was not found')
|
---|
83 | self.f = ROOT.factfits(path)
|
---|
84 | self._make_header()
|
---|
85 | self._setup_columns()
|
---|
86 |
|
---|
87 | def _make_header(self):
|
---|
88 | """
|
---|
89 | """
|
---|
90 | str_to_bool = { 'T':True, 'F':False}
|
---|
91 | type_conversion = { 'I' : int, 'F' : float, 'T' : str, 'B' : str_to_bool.__getitem__}
|
---|
92 |
|
---|
93 | self.header = {}
|
---|
94 | self.header_comments = {}
|
---|
95 | for key,entry in self.f.GetKeys():
|
---|
96 | try:
|
---|
97 | self.header[key] = type_conversion[entry.type](entry.value)
|
---|
98 | except KeyError:
|
---|
99 | raise IOError("Error: entry type unknown.\n Is %s, but should be one of: [I,F,T,B]" % (entry.type) )
|
---|
100 | self.header_comments[key] = entry.comment
|
---|
101 |
|
---|
102 | def _setup_columns(self):
|
---|
103 | """
|
---|
104 | """
|
---|
105 | col_type_to_np_type_map = { 'L' : 'b1', 'A' : 'a1', 'B' : 'i1',
|
---|
106 | 'I' : 'i2', 'J' : 'i4', 'K' : 'i8', 'E' : 'f4', 'D' : 'f8'}
|
---|
107 | self.cols = {}
|
---|
108 | for key,col in self.f.GetColumns():
|
---|
109 | self.cols[key] = np.zeros(col.num, col_type_to_np_type_map[col.type])
|
---|
110 | if col.num != 0:
|
---|
111 | self.f.SetPtrAddress(key, self.cols[key])
|
---|
112 |
|
---|
113 | def __iter__(self):
|
---|
114 | return self
|
---|
115 |
|
---|
116 | def next(self, row=None):
|
---|
117 | """
|
---|
118 | """
|
---|
119 | if row is None:
|
---|
120 | if self.f.GetNextRow() == False:
|
---|
121 | raise StopIteration
|
---|
122 | else:
|
---|
123 | row = int(row)
|
---|
124 | if self.f.GetRow(row) == False:
|
---|
125 | raise StopIteration
|
---|
126 | return self
|
---|
127 |
|
---|
128 | class AuxFile( Fits ):
|
---|
129 | """ easy(?) access to FACT aux files
|
---|
130 | """
|
---|
131 | __module__ = 'pyfact'
|
---|
132 | def __init__(self, path, verbose=False):
|
---|
133 | self._verbose = verbose
|
---|
134 | super(AuxFile, self).__init__(path)
|
---|
135 |
|
---|
136 | def _setup_columns(self):
|
---|
137 | col_type_to_np_type_map = { 'L' : 'b1', 'A' : 'a1', 'B' : 'i1',
|
---|
138 | 'I' : 'i2', 'J' : 'i4', 'K' : 'i8', 'E' : 'f4', 'D' : 'f8'}
|
---|
139 | self._cols = {}
|
---|
140 | self.cols = {}
|
---|
141 | N = self.header['NAXIS2']
|
---|
142 | for key,col in self.f.GetColumns():
|
---|
143 | self._cols[key] = np.zeros(col.num, col_type_to_np_type_map[col.type])
|
---|
144 | self.cols[key] = np.zeros((N, col.num), col_type_to_np_type_map[col.type])
|
---|
145 | if col.num != 0:
|
---|
146 | self.f.SetPtrAddress(key, self._cols[key])
|
---|
147 |
|
---|
148 |
|
---|
149 | for i,row in enumerate(self):
|
---|
150 | if self._verbose:
|
---|
151 | try:
|
---|
152 | step = int(self._verbose)
|
---|
153 | except:
|
---|
154 | step = 10
|
---|
155 | if i % step == 0:
|
---|
156 | print "reading line", i
|
---|
157 |
|
---|
158 | for key in self._cols:
|
---|
159 | self.cols[key][i,:] = self._cols[key]
|
---|
160 |
|
---|
161 |
|
---|
162 |
|
---|
163 |
|
---|
164 | class RawData( Fits ):
|
---|
165 | """ Special raw data FITS file access (with DRS4 calibration)
|
---|
166 |
|
---|
167 | During iteration the C++ method DrsCalibration::Apply is being called.
|
---|
168 | """
|
---|
169 | __module__='pyfact'
|
---|
170 | def __init__(self, data_path, calib_path):
|
---|
171 | """ -constructor-
|
---|
172 | *data_path* : fits or fits.gz file of the data including the path
|
---|
173 | *calib_path* : fits or fits.gz file containing DRS calibration data
|
---|
174 | """
|
---|
175 | super(RawData, self).__init__(data_path)
|
---|
176 | self.cols['CalibData'] = np.zeros( self.cols['Data'].shape, np.float32)
|
---|
177 | self.cols['CalibData2D'] = self.cols['CalibData'].reshape( self.header['NPIX'], -1)
|
---|
178 | if not self.cols['CalibData2D'].base is self.cols['CalibData']:
|
---|
179 | print "Error seomthing went wrong!"
|
---|
180 | self.drs_calibration = ROOT.DrsCalibration()
|
---|
181 | self.drs_calibration.ReadFitsImp( calib_path )
|
---|
182 |
|
---|
183 | self.drs_calibrate = ROOT.DrsCalibrate()
|
---|
184 | self.list_of_previous_start_cells = []
|
---|
185 |
|
---|
186 | def __iter__(self):
|
---|
187 | """ iterator """
|
---|
188 | return self
|
---|
189 |
|
---|
190 | def next(self, row=None):
|
---|
191 | """
|
---|
192 | """
|
---|
193 | super(RawData, self).next(row)
|
---|
194 |
|
---|
195 | self.drs_calibration.Apply( self.cols['CalibData'],
|
---|
196 | self.cols['Data'],
|
---|
197 | self.cols['StartCellData'],
|
---|
198 | self.header['NROI'])
|
---|
199 |
|
---|
200 | for previous_start_cells in self.list_of_previous_start_cells:
|
---|
201 | self.drs_calibrate.CorrectStep(
|
---|
202 | self.cols['CalibData'],
|
---|
203 | self.header['NPIX'],
|
---|
204 | self.header['NROI'],
|
---|
205 | previous_start_cells,
|
---|
206 | self.cols['StartCellData'],
|
---|
207 | self.header['NROI']+10)
|
---|
208 | self.drs_calibrate.CorrectStep(
|
---|
209 | self.cols['CalibData'],
|
---|
210 | self.header['NPIX'],
|
---|
211 | self.header['NROI'],
|
---|
212 | previous_start_cells,
|
---|
213 | self.cols['StartCellData'],
|
---|
214 | 3)
|
---|
215 | self.list_of_previous_start_cells.append(self.cols['StartCellData'])
|
---|
216 | if len(self.list_of_previous_start_cells) > 5:
|
---|
217 | self.list_of_previous_start_cells.pop(0)
|
---|
218 |
|
---|
219 | for ch in range(self.header['NPIX']):
|
---|
220 | self.drs_calibrate.RemoveSpikes3(self.cols['CalibData2D'][ch], self.header['NROI'])
|
---|
221 |
|
---|
222 | return self
|
---|
223 |
|
---|
224 |
|
---|
225 |
|
---|
226 | if __name__ == '__main__':
|
---|
227 | """ Example """
|
---|
228 | print "Example for calibrated raw-file"
|
---|
229 | f = RawData(sys.argv[1], sys.argv[2])
|
---|
230 |
|
---|
231 | print "number of events:", f.header['NAXIS2']
|
---|
232 | print "date of observation:", f.header['DATE']
|
---|
233 | print "The files has these cols:", f.cols.keys()
|
---|
234 |
|
---|
235 | for counter,row in enumerate(f):
|
---|
236 | print "Event Id:", row.cols['EventNum']
|
---|
237 | print "shape of column 'StartCellData'", row.cols['StartCellData'].shape
|
---|
238 | print "dtype of column 'Data'", row.cols['StartCellData'].dtype
|
---|
239 | if counter > 3:
|
---|
240 | break
|
---|
241 | # get next row
|
---|
242 | f.next()
|
---|
243 | print "Event Id:", f.cols['EventNum']
|
---|
244 | # get another row
|
---|
245 | f.next(10)
|
---|
246 | print "Event Id:", f.cols['EventNum']
|
---|
247 | # Go back again
|
---|
248 | f.next(3)
|
---|
249 | print "Event Id:", f.cols['EventNum']
|
---|
250 |
|
---|
251 | import matplotlib.pyplot as plt
|
---|
252 | plt.ion()
|
---|
253 | """
|
---|
254 | for i in range(f.header['NPIX']):
|
---|
255 | plt.cla()
|
---|
256 | plt.title("Event %d"%(f.cols['EventNum']))
|
---|
257 | plt.plot(f.cols['CalibData2D'][i], '.:', label='pixel %d'%i)
|
---|
258 | plt.legend()
|
---|
259 | answer = raw_input('anykey for next pixel; "q" to quit. :')
|
---|
260 | if 'q' in answer.lower():
|
---|
261 | break
|
---|
262 | """
|
---|
263 |
|
---|
264 |
|
---|
265 |
|
---|
266 | from scipy import signal
|
---|
267 | from scipy.signal import freqs, iirfilter
|
---|
268 | #b, a = iirfilter(4, [1, 10], 1, 60, analog=True, ftype='cheby1')
|
---|
269 | #b, a = signal.iirfilter(17, [50, 200], rs=60, btype='band',
|
---|
270 | # analog=True, ftype='cheby2')
|
---|
271 | #w, h = signal.freqs(b, a, 1000)
|
---|
272 | #w, h = freqs(b, a, worN=np.logspace(-1, 2, 1000))
|
---|
273 |
|
---|
274 | b, a = iirfilter(N=4, # the order of the filter
|
---|
275 | Wn = [0.15, 0.25], # critical frequencies
|
---|
276 | rp = 1, # maximum ripple in passband
|
---|
277 | rs = 60, # minimum attenuation in stopband
|
---|
278 | btype='bandstop', # {'bandpass', 'lowpass', 'highpass', 'bandstop'}
|
---|
279 | #analog=True,
|
---|
280 | ftype='cheby1')
|
---|
281 | #w, h = freqs(b, a, worN=np.logspace(-1, 5, 1000))
|
---|
282 | w, h = freqs(b, a)
|
---|
283 | plt.figure()
|
---|
284 | plt.semilogx(w, 20 * np.log10(abs(h)))
|
---|
285 | plt.xlabel('Frequency')
|
---|
286 | plt.ylabel('Amplitude response [dB]')
|
---|
287 | plt.grid()
|
---|
288 | plt.show()
|
---|
289 |
|
---|
290 |
|
---|
291 |
|
---|
292 |
|
---|
293 |
|
---|
294 | fig1, (ax11,ax12) = plt.subplots(nrows=2)
|
---|
295 | fig2, (ax21,ax22) = plt.subplots(nrows=2)
|
---|
296 |
|
---|
297 | flt = firwin(13, [0.15,0.25])
|
---|
298 | t = np.arange(f.cols['CalibData2D'].shape[1])/2.*1e-9
|
---|
299 | fig3, (ax31, ax32) = plt.subplots(nrows=2)
|
---|
300 | ax31.plot(flt)
|
---|
301 | #w, h = freqs(flt, [1.0])
|
---|
302 | w, h = freqs(flt, [1.0], worN=np.logspace(-1, 5, 1000))
|
---|
303 | ax32.semilogx(w, 20 * np.log10(abs(h)))
|
---|
304 | ax32.set_xlabel('Frequency')
|
---|
305 | ax32.set_ylabel('Amplitude response [dB]')
|
---|
306 | ax32.grid()
|
---|
307 |
|
---|
308 |
|
---|
309 |
|
---|
310 | for i in range(f.header['NPIX']):
|
---|
311 | data = f.cols['CalibData2D'][i]
|
---|
312 | fil_data = lfilter(flt, [1.0], data)
|
---|
313 |
|
---|
314 | ax11.cla()
|
---|
315 | ax12.cla()
|
---|
316 | ax11.set_title("Event %d"%(f.cols['EventNum']))
|
---|
317 | ax11.plot(t, data, '.:', label='pixel %d'%i)
|
---|
318 | ax11.legend()
|
---|
319 | ax11.set_ylim(-20,40)
|
---|
320 | ax12.plot(t[:fil_data.size], fil_data, '.:', label='pixel %d'%i)
|
---|
321 | ax12.legend()
|
---|
322 | ax12.set_ylim(-20,40)
|
---|
323 |
|
---|
324 | """
|
---|
325 | Pxx, freqs, t, plot = pylab.specgram( f.cols['CalibData2D'][i],
|
---|
326 | NFFT=32, Fs=2000000000,
|
---|
327 | detrend=pylab.detrend_none,
|
---|
328 | window=pylab.window_hanning,
|
---|
329 | noverlap=16)
|
---|
330 | """
|
---|
331 | ax21.cla()
|
---|
332 | ax22.cla()
|
---|
333 | dat, freqs, bins, im = ax21.specgram(data,
|
---|
334 | NFFT=64, Fs=2000000000,
|
---|
335 | detrend=pylab.detrend_none,
|
---|
336 | window=pylab.window_hanning,
|
---|
337 | noverlap=60,
|
---|
338 | interpolation='nearest')
|
---|
339 | ax21.axis('tight')
|
---|
340 | ax21.set_ylabel("frequency [Hz]")
|
---|
341 | ax21.set_xlabel("time [s]")
|
---|
342 |
|
---|
343 | dat2, freqs2, bins2, im2 = ax22.specgram(data,
|
---|
344 | NFFT=64, Fs=2000000000,
|
---|
345 | detrend=pylab.detrend_none,
|
---|
346 | window=pylab.window_hanning,
|
---|
347 | noverlap=60,
|
---|
348 | interpolation='nearest')
|
---|
349 | ax22.axis('tight')
|
---|
350 | ax22.set_ylabel("frequency [Hz]")
|
---|
351 | ax22.set_xlabel("time [s]")
|
---|
352 |
|
---|
353 | plt.show()
|
---|
354 |
|
---|
355 | answer = raw_input('anykey for next pixel; "q" to quit. :')
|
---|
356 | if 'q' in answer.lower():
|
---|
357 | break
|
---|